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Abstract 

The  shallow  water  acoustic  channel  supports  far-field  propagation  in  a  discrete  set  of 
modes.  Ocean  experiments  have  confirmed  the  modal  nature  of  acoustic  propagation, 
but  no  experiment  has  successfully  excited  only  one  of  the  suite  of  mid-frequency 
propagating  modes  propagating  in  a  coastal  environment.  The  ability  to  excite  a 
single  mode  would  be  a  powerful  tool  for  investigating  shallow  water  ocean  pro¬ 
cesses.  A  feedback  control  algorithm  incorporating  elements  of  adaptive  estimation, 
underwater  acoustics,  array  processing  and  control  theory  to  generate  a  high-fidelity 
single  mode  is  presented.  This  approach  also  yields  a  cohesive  framework  for  evalu¬ 
ating  the  feasibility  of  generating  a  single  mode  with  given  array  geometries,  noise 
characteristics  and  source  power  limitations.  Simulations  and  laboratory  waveguide 
experiments  indicate  the  proposed  algorithm  holds  promise  for  ocean  experiments. 
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Chapter  1 


Introduction 


The  shallow  water  channel  is  a  waveguide  for  pressure  waves,  supporting  propagation 
for  a  discrete  set  of  normal  modes  in  the  far  field,  [1],  [2].  Ocean  acoustics  problems 
such  as  internal  wave  tomography  [3],  [4]  remote  sensing  [5],  and  noise  propagation 
[6]  [7]  often  discuss  propagation  in  terms  of  modes.  Numerous  experiments  have 
verified  the  modal  model  of  propagation  in  the  shallow  water  [8],  [9],  [2],  [10],  [11], 
but  no  mid-frequency  (circa  400  Hz)  ocean  experiment  has  successfully  excited  a 
single  mode  in  a  shallow  water  environment. 

This  thesis  proposes  an  algorithm  for  controlling  a  vertical  array  of  narrowband 
sources  such  that  the  pressure  field  in  the  shallow  water  acoustic  channel  measured 
at  a  reference  location  consists  of  only  a  single  mode.  While  none  of  the  individ¬ 
ual  aspects  of  the  proposed  algorithm  are  original,  the  synthesis  of  these  elements 
produces  a  novel  approach  to  exciting  a  single  mode.  The  algorithm  uses  feedback 
control  to  obtain  the  desired  pressure  field  at  the  reference  hydrophone  array.  This 
technique  is  commonly  used  for  pressure  field  control  in  open-air  acoustics  applica¬ 
tions  such  as  active  noise  control  [12],  [13]  but  is  not  often  used  in  the  underwater 
acoustics  community.  The  specific  feedback  technique  used  to  control  the  pressure 
field  at  the  feedback  array  is  the  method  of  indirect  control  as  described  by  Narendra 
and  Annaswamy  [14].  The  control  algorithm  estimates  the  Green’s  function  between 
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each  element  of  the  source  array  and  each  hydrophone  of  the  reference  array.  This 
matrix  of  estimated  Green’s  functions  is  then  inverted  to  determine  the  shading  to 
apply  to  the  source  array  in  order  to  obtain  the  desired  pressure  field  at  the  feedback 
hydrophone  array.  Both  the  acoustic  model  and  control  algorithm  are  simple  and 
standard  techniques  in  their  respective  fields.  Although  much  more  sophisticated 
approaches  to  the  feedback  control  algorithm  and  ocean  acoustics  are  possible,  this 
simple  approach  will  allow  a  clearer  understanding  of  the  problem  in  this  preliminary 
investigation.  If  the  algorithm  proposed  using  these  straightforward  ideas  proves  suc¬ 
cessful  here,  further  investigations  can  determine  if  more  sophisticated  approaches 
to  both  the  acoustics  and  control  can  provide  further  gains. 

Although  the  feedback  control  algorithm  requires  measurements  of  the  sound 
speed  profile  at  the  feedback  array,  it  obviates  the  need  for  detailed  a  priori  infor¬ 
mation  throughout  the  control  volume  in  order  to  excite  the  desired  pressure  field. 
Previous  work  exciting  a  single  mode  in  a  laboratory  waveguide  exploited  detailed 
environmental  knowledge  and  time-invariance  to  find  the  array  weights  for  exciting 
a  single  mode  analytically  [15],  [16].  While  these  experiments  successfully  controlled 
the  pressure  in  laboratory  tanks  in  open  loop  mode,  this  approach  appears  unlikely 
to  work  in  an  arbitrary  ocean  situation  without  detailed  knowledge  of  the  environ¬ 
ment.  Range  inhomogeneities  in  the  environment  may  couple  energy  among  modes 
or  from  the  continuum  into  propagating  modes  before  the  pressure  wave  reaches  the 
desired  observation  volume. 

The  distribution  of  propagating  modes  can  be  inferred  from  the  pressure  field 
sampled  at  the  locations  of  the  hydrophones  in  the  reference  array.  This  estimation 
problem  is  known  as  mode  filtering  in  the  ocean  acoustics  literature  [17],  [9].  A 
variety  of  algorithms  have  been  proposed  for  solving  this  inverse  problem.  The 
sensitivity  of  this  estimate  to  noise  depends  on  the  geometry  of  the  hydrophone  array 
and  the  mode  filtering  algorithm.  A  crucial  issue  for  this  thesis  is  determining  how 
the  error  between  the  desired  and  observed  pressure  profiles  is  related  to  the  error 
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between  the  desired  and  actual  mode  distribution  propagating  in  the  channel.  This 
problem  and  mode  filtering  are  intrinsically  related.  The  control  theory  framework 
developed  in  this  thesis  for  determining  the  observability  of  modes  can  offer  insight 
into  this  problem. 


Figure  1-1:  Proposed  Experimental  Setup 

Figure  1-1  shows  the  proposed  experiment  configuration.  The  region  between 
the  source  and  feedback  arrays  will  be  called  the  feedback  volume,  while  the  region 
between  the  feedback  and  observation  arrays  will  be  referred  to  as  the  observation 
volume.  The  feedback  algorithm  attempts  to  control  the  pressure  field  such  that 
only  a  single  mode  is  propagating  at  the  start  of  the  observation  volume,  i.e.,  at  the 
feedback  array.  This  allows  the  observers  at  the  observation  array  to  be  confident 
the  pressure  field  they  measure  was  generated  by  a  high-fidelity  mode  impinging  on 
the  observation  volume.  The  source  array  is  tended  by  a  ship,  while  the  feedback 
array  is  located  at  rp  =  1-2  km  downrange,  at  the  start  of  the  far  field  [18].  The 
far-field  is  defined  as  beginning  at  the  range  where  all  the  significant  energy  prop¬ 
agating  in  the  waveguide  is  described  by  the  trapped  modes.  The  feedback  array 
transmits  the  pressure  observed  at  its  hydrophones  back  to  the  source  ship  over  a  ra- 
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dio  telemetry  link.  The  scientific  observation  array  is  deployed  10-20  km  downrange 
to  measure  the  distribution  of  modes  emerging  from  the  observation  region.  For  the 
purposes  of  the  pressure  received  at  the  observation  array,  the  feedback  array  can  be 
considered  a  virtual  source  exciting  the  ocean  at  that  location  with  a  single  mode. 
Knowing  the  pressure  field  entering  the  observation  volume  consisted  of  only  a  single 
mode,  the  modes  arriving  at  the  observation  array  can  give  information  about  the 
oceanographic  properties  of  the  intervening  water  mass. 

Several  studies  have  examined  the  effect  of  internal  waves  on  acoustic  propagation 
in  the  coastal  environment.  Lee  [19]  used  a  ray  approach  to  demonstrate  the  presence 
of  internal  waves  supported  on  a  thermocline  intensified  the  contrast  in  the  shadow 
zones  compared  with  a  similar  environment  without  the  internal  waves.  Studies  by 
Zhou,  Zhang  and  Rogers  [4]  and  Peregrym  [3]  simulated  the  propagation  of  acoustic 
normal  modes  through  internal  wave  packets  with  sinusoidal  displacement  profiles. 
These  simulations  demonstrated  a  narrowband  resonance  coupling  energy  between 
modes  when  the  difference  between  the  modes’  horizontal  wavenumbers  equaled  the 
spatial  wavenumber  of  the  internal  wave  packet.  More  recent  work  by  Preisig  and 
Duda  [20]  have  found  the  earlier  studies  may  have  exaggerated  the  narrowness  of 
the  resonance  by  using  an  unrealistically  periodic  displacement  profile.  All  of  this 
work  indicates  that  given  the  pressure  field  entering  the  observation  volume  contains 
only  one  mode,  the  modes  observed  at  the  observation  array  contain  significant 
information  about  the  properties  of  any  internal  waves  in  that  volume. 

The  problem  of  exciting  a  single  mode  using  feedback  control  incorporates  ele¬ 
ments  of  underwater  acoustics,  mode  filtering,  control  theory,  and  signal  processing. 
The  subsequent  sections  of  this  chapter  summarize  relevant  aspects  of  each  of  these 
disciplines. 

Chapter  2  describes  the  proposed  control  algorithm  obtained  by  synthesizing  the 
material  covered  in  Chapter  1.  First,  a  model  is  proposed  for  the  acoustic  channel 
of  the  feedback  volume.  The  indirect  control  algorithm  developed  in  this  chapter 
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first  estimates  the  parameters  of  the  channel  model,  then  uses  this  model  to  choose 
the  complex  source  weights  to  excite  the  desired  pressure  field  at  the  feedback  array. 
Several  different  estimators  are  presented  for  the  identification  of  the  channel  model. 

The  proposed  control  algorithm  is  evaluated  in  Chapter  3  in  simulations  of  a  vari¬ 
ety  of  shallow  water  ocean  environments.  All  of  the  environments  modeled  are  based 
on  measured  profiles  from  the  South  Continental  Shelf  off  Martha’s  Vineyard.  The 
propagation  for  these  environments  is  simulated  using  the  finite-element  parabolic 
equation  (FEPE)  approximation  to  the  wave  equation  [21]. 

In  addition  to  the  simulations,  the  thesis  presents  results  from  a  series  of  ex¬ 
periments  using  the  algorithm  to  control  the  pressure  field  at  a  hydrophone  array 
in  a  scale  model  laboratory  waveguide.  These  experiments,  described  in  Chapter 
4,  demonstrate  the  algorithm  is  robust  enough  to  work  in  real  time  with  actual 
sources,  hydrophones,  and  acoustic  waves  propagating  through  water.  The  different 
algorithms  for  estimating  the  channel  model  are  examined  for  both  steady  state  and 
transient  performance. 

Chapter  5  makes  conclusions  about  the  algorithm  proposed  based  on  the  simu¬ 
lations  and  experiments  described  in  the  thesis.  In  addition,  this  chapter  suggests 
future  directions  for  developing  algorithms  for  single  mode  excitation  and  other  pos¬ 
sible  uses  for  feedback  control  and  single  mode  excitation  in  ocean  acoustics. 


1.1  Normal  Mode  Model  for  Acoustic  Propaga¬ 
tion 


The  pressure  field  for  a  time-harmonic  point  source  located  at  fq  in  a  range-independent 
environment  can  be  described  by  the  wave  equation  [22], 


p(z)V 


1 

p{z) 


Vp(r) 


+  k^{z)p{r:)  =  -47r5(r  -  fq). 


(1.1) 
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where  r  =  (r,  6,  z)  is  the  observer’s  location  in  cylindrical  coordinates,  Tq  =  (0, 0, 2:0) 
is  the  source  location,  and  the  wavenumber  k{z)  =  ujfc^z)  is  the  ratio  of  the  angular 
frequency  to  the  depth-dependent  sound  speed.  Assuming  cylindrically  symmetric 
solutions  separable  in  range  and  depth  of  the  form 

p(r,z)  =  ^',n(2;)i4i(r)  (1.2) 


yields  the  normal  mode  solution  for  the  field  [22],  [23],  and  [1].  Substituting  Eq.  1.2 
into  the  homogeneous  form  of  Eq.  1.1  and  separating  the  functions  of  depth  and 
range  gives  the  differential  equations  specifying  the  eigenfunctions  for  each  of  these 
coordinates: 


1 

p{z)  dz^ 


+ 


dz  [p{z) 
r  dr 


+ 


,  klJz) 

dz 

dRmir) 


dr 


P(^) 


0 

0, 


(1.3) 

(1.4) 


where  the  separation  constant  is  k^,  its  square  root  km  is  the  horizontal  wavenumber, 
and  k^m{z)  =  y/k^^iz)  -  is  the  vertical  wavenumber. 

The  total  solution  for  the  pressure  will  be  a  superposition  of  all  solutions  of  the 
form  of  Eq.  1.2 

p(^  0'm{Zo)'^m{z)Rm{r)-  (1-5) 

m 

Substituting  this  expression  into  Eq.  1.1  and  simplifying  the  resulting  equation  with 
Eq.  1.3  reveals 


Rm{r)  =  iTTH^^\kmr)  (1.6) 

am{Zo)  =  '^m{zo)/p{zQ),  (1.7) 

where  H^\kmT)  is  the  zeroth-order  Hankel  function  of  the  first  kind,  and  the  oper¬ 
ator  (•)*  denotes  complex  conjugation.  Combining  these  gives  the  solution  for  the 
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pressure  field: 


(1.8) 


p{r,  z)  =  ^*mi^)^m{z)H^^\kmr). 

For  the  far  field  {kmr  »  1),  the  Hankel  function  can  be  replaced  by  its  asymptotic 
approximation  to  get 

Jknir 

Y.^m^^o)^rn{z)^=.  (1.9) 

P('^o)  m  ykfnX 

In  general,  the  sound  speed  profile  c  and  density  p  are  functions  of  range  as  well 
as  depth.  Consequently,  the  solutions  to  Eq.  1.3  vary  with  range.  The  solutions  to 
the  depth  eigenfunction  equation  for  the  hydrographic  and  boundary  conditions  at  a 
given  range  are  known  as  the  local  modes  at  that  range.  The  pressure  field  at  range 
r  can  be  described  as  a  superposition  of  the  local  modes  "^miz]  r),  or 

p{r,  z)  =  Y^  dm{r,  zo)'^miz-,  r),  (1.10) 

m 

where  dm(r,Zo)  are  the  complex  mode  coefficients.  If  the  channel  is  linear,  the 
pressure  field  excited  by  an  array  of  sources  can  be  written  as  a  superposition  of 
solutions  of  the  form  of  Eq.  1.10.  Equivalently,  the  complex  mode  coefficients  can 
be  considered  to  be  functions  of  the  source  depth  vector  Zg  and  the  complex  source 

weights  w,  in  addition  to  r.  The  pressure  field  generated  by  a  vertical  array  of 

sources  at  depths  z*  weighted  by  w  is 

P(r,  z)  =  J2  ^rn{r,  z„  w)^,„(2;;  r).  (1.11) 

m 

The  proposed  algorithm  attempts  to  choose  w  such  that  only  one  mode  coefficient 
dmo  is  non-zero  at  the  feedback  hydrophone  array  at  r  =  r^. 

Acoustic  propagation  of  range-dependent  modes  is  categorized  as  adiabatic  or 
coupled.  The  adiabatic  model  of  propagation  assumes  no  energy  is  exchanged  be¬ 
tween  different  modes  although  the  local  modes  '^m{z',  r)  may  vary  slowly  with  range. 
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In  coupled  propagation,  energy  is  exchanged  between  the  local  modes  due  to  range 
variations  in  bathymetry  and  sound  speed  [24], [25]. 

The  normal  mode  formulation  for  acoustic  propagation  is  especially  attractive 
in  shallow  water  environments  at  medium  to  low  frequency.  In  this  scenario,  the 
horizontal  wavenumber  spectrum  can  be  segmented  into  three  classes  of  waves.  The 
propagating,  or  “trapped,”  modes  are  the  solutions  to  the  eigenfunction  equations 
(Eqs.  1.3  and  1.4)  such  that  the  horizontal  wavenumber  km  is  real,  or  at  least  has  only 
a  very  small  imaginary  part  due  to  medium  attenuation  during  propagation.  The 
evanescent  modes  have  vertical  profiles  which  are  solutions  to  the  depth  eigenfunction 
equation  (Eq.  1.3),  but  with  k;sm  >  k,  so  km  =  -  kzm  is  imaginary.  These 

modes  decay  exponentially  in  range.  In  an  ideal  waveguide  with  perfectly  reflecting 
boundaries,  the  trapped  and  evanescent  modes  contain  all  the  energy.  For  a  realistic 
ocean  waveguide,  the  situation  is  complicated  by  the  absence  of  a  perfectly  reflecting 
bottom.  The  solutions  in  this  environment  can  be  found  by  taking  the  limit  as  the 
ideal  boundary  goes  to  infinite  depth  [22].  This  results  in  a  continuum  of  horizontal 
wavenumbers  from  the  solutions  to  Eqs.  1.3  and  1.4.  A  discrete  finite  subset  of  these 
solutions  propagate  in  range  with  little  or  no  attenuation;  these  solutions  are  the 
trapped  modes  for  the  waveguide.  The  rest  of  the  solutions  constitute  the  modal 
continuum.  The  contribution  of  the  continuum  is  significant  at  short  ranges,  but 
quickly  decays  in  range  to  leave  only  the  trapped  modes  in  the  far  field. 

Many  shallow  water  environments  have  bathymetry  that  is  sufficiently  range- 
varying  to  invalidate  the  adiabatic  propagation  model.  Desaubies  [25]  quantified 
measures  indicating  the  importance  of  modal  coupling  in  an  environment.  He  found 
that  the  precise  conditions  when  the  adiabatic  approximation  breaks  down  depends 

on  the  quantity  of  interest,  i.e.,  pressure,  transmission  loss,  propagation  time,  etc _ 

For  many  of  the  quantities  of  interest,  Desaubies  found  the  conditions  when  the 
adiabatic  approximation  became  invalid  depended  in  subtle  and  sometimes  intricate 
ways  on  the  rate  of  change  of  the  environment  in  range.  Because  the  coupling  induced 
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by  the  bottom  or  range-inhomogeneities  in  the  water  column  can  be  strong,  it  is  very 
difficult  to  compute  the  source  weights  w  to  excite  only  a  single  mode  without  very 
extensive  and  accurate  environmental  measurements. 

Many  researchers  have  investigated  and  validated  the  modal  model  of  shallow 
water  propagation.  The  seminal  work  in  shallow  water  modal  propagation  is  the 
monograph  “Theory  of  Propagation  of  Explosive  Sound  in  Shallow  Water”  by  C.  L. 
Pekeris  [2].  This  paper  develops  the  so-called  “Pekeris  waveguide,”  consisting  of  an 
isovelocity  water  layer,  bounded  above  by  a  pressure  release  boundary  and  below  by 
a  higher-velocity  isovelocity  halfspace.  Using  this  model  and  the  theory  of  normal 
modes,  the  paper  predicts  the  features  of  the  pressure  pulse  caused  by  an  impulsive 
source  (TNT  detonation)  on  the  ocean  bottom  [26]. 

Bucker  did  significant  work  verifying  the  modal  model  of  propagation  for  the 
shallow  water  channel.  In  [10],  he  compared  experimental  and  theoretical  mode 
curves  for  explosive  sources  in  the  Bering  Sea,  obtaining  close  agreement.  The  time- 
frequency  distribution  of  the  arrivals  allowed  him  to  identify  the  various  modes 
in  the  energy  observed  at  a  single  hydrophone.  In  a  later  paper  [27],  he  computed 
expressions  for  the  normal  mode  shapes  for  a  variety  of  analytic  sound  speed  profiles, 
as  well  as  one  measured  experimental  environment.  Bucker  used  the  computed  mode 
shapes  to  predict  transmission  loss  as  a  function  of  range  in  that  experiment,  and 
found  reasonable  agreement  with  the  observed  values. 

Ferris  et  al.  [11],  [8],  transmitted  time- windowed  sinusoidal  pulses  from  a  single 
400  Hz  source  to  a  nine  element  vertical  array  of  hydrophones  10  km  downrange. 
The  differing  group  delays  of  the  modes  resulted  in  different  travel  times  between  the 
source  and  receiver.  These  travel  times  differed  by  enough  that  two,  and  sometimes 
three,  distinct  modes  could  be  temporally  resolved  at  the  receiver  without  any  spatial 
processing.  The  vertical  profiles  of  these  modes  across  the  array  corresponded  well 
with  the  computed  mode  shapes  for  the  sound  speed  profile.  The  later  paper,  [8],  also 
employed  a  simple  spatial  filter  which  projected  the  pressure  field  against  samples 
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of  the  theoretically  derived  mode  shapes  to  verify  the  mode  strength  as  a  function 
of  source  depth  matched  the  computed  mode  shapes. 

Tindle  et  al.  [9]  performed  an  experiment  very  clearly  demonstrating  shallow 
water  propagation  was  well-modeled  by  modes.  They  found  a  site  almost  perfectly 
homogeneous  in  range  and  depth.  Over  a  range  of  5  km,  the  bottom  depth  was 
50  ±  1  m,  with  a  sound  speed  profile  that  was  nearly  perfectly  homogeneous  in  depth 
and  range  with  value  1508.7±0.3  m/s.  The  analytically  derived  depth  eigenfunctions 
for  an  isovelocity  water  column  fit  this  experimental  environment  very  well.  By 
running  experiments  transmitting  windowed  sinusoids  centered  at  60,  100,  and  140 
Hz  and  using  the  least  squares  mode  filter  which  is  discussed  in  Section  1.3,  they 
were  able  to  resolve  the  predicted  arrivals  of  one,  two  and  three  modes,  respectively. 

Buckingham  [28],  [29]  analytically  derived  an  approximate  expression  for  the 
modes  for  a  downslope  wedge  environment  with  an  isovelocity  water  column  over  an 
isovelocity  halfspace  bottom.  Although  the  propagation  for  this  environment  can  be 
approximated  with  the  coupled  mode  model,  he  demonstrated  the  true  eigenfunctions 
of  the  wave  equation  for  this  geometry  were  circularly  curved  wavefronts,  centered  on 
the  apex  of  the  wedge,  sinusoidally  varying  in  angle.  These  wedge  modes  propagate 
downslope  without  coupling  energy  among  them,  and  thus  are  the  true  modes  of 
the  system.  Since  the  wedge  possesses  spherical  symmetry  rather  cylindrical,  it 
is  intuitively  reasonable  that  the  azimuthal  angle  plays  the  role  depth  played  in 
cylindrical  geometry  with  a  flat  bottom. 

Tindle,  Hobaek,  and  Muir  [30],  [31]  confirmed  Buckingham’s  theory  in  a  series  of 
thorough  and  elegant  experiments  investigating  modal  propagation  in  a  scale  model 
downslope  environment.  Their  laboratory  tank  was  93  cm  wide  and  10  m  long  with 
10  cm  of  water  at  the  source  over  a  20  cm  deep  sand  bottom.  The  bottom  could 
be  tilted  at  any  angle  between  0°  and  9°.  Their  source  transmitted  pulses  centered 
on  80  kHz  and  the  pressure  field  in  the  tank  was  measured  by  a  hydrophone  moved 
between  successive  transmissions  to  create  a  synthetic  discrete  array.  Tindle  et  al. 
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observed  the  signals  propagated  downslope  with  curved  wavefronts  centered  on  the 
apex  of  the  wedge  as  predicted  by  Buckingham.  They  resolved  these  modes  by 
spatially  filtering  the  pressure  field  using  the  least-squares  technique  of  Tindle  et 
al.  in  [9].  Comparing  the  results  obtained  using  a  mode  filter  designed  with  the 
local  vertical  modes  to  those  obtained  when  the  filter  used  the  curved  wedge  modes 
demonstrated  that  only  the  latter  propagated  without  coupling. 

Several  experiments  examined  single  mode  excitation  in  laboratory  tanks.  Clay 
and  Huang  attempted  to  transmit  and  receive  a  single  mode  in  a  laboratory  tank  as 
part  of  an  experiment  to  measure  backscattering  from  fish  [15].  The  tank  they  used 
had  25  cm  of  water  and  they  transmitted  at  a  frequency  of  220  Hz.  The  experimen¬ 
tal  waveguide  had  pressure-release  surfaces  both  above  and  below,  thus  supported 
roughly  70  modes  at  this  frequency.  Clay  and  Huang  attempted  to  excite  only  the 
first  mode  using  8  source  transducers,  shading  these  sources  with  samples  of  the 
desired  mode  shape  quantized  to  three  levels.  The  same  array  and  shading  were 
used  as  a  receiver.  Under  these  conditions  of  gross  undersampling  and  coarse  quan¬ 
tization,  it  is  not  surprising  that  the  received  pressure  field  roughly  matched  the 
shading.  While  Clay  and  Huang  claim  this  indicates  only  the  first  mode  was  excited, 
it  also  seems  possible  they  synthesized  an  approximation  to  this  crudely  quantized 
first  mode  shape  using  all  the  modes  as  basis  functions.  Thus,  the  channel  may 
not  have  contained  just  one  mode,  but  may  actually  have  contained  many  propagat¬ 
ing  modes  whose  superposition  approximated  a  quantized  version  of  the  first  mode 
shape.  Because  Clay  and  Huang’s  receiver  undersampled  the  channel  spatially  and 
they  gave  no  information  about  temporal  dispersion,  it  is  not  possible  to  conclude 
with  certainty  from  their  results  what  modes  were  propagating  in  the  channel.  Also, 
their  experiment  used  fixed  array  weights  built  into  the  source/receiver  array.  The 
laboratory  environment  was  almost  certainly  time-invariant,  so  the  fixed  nature  of 
the  weights  was  probably  not  an  issue.  Because  it  could  not  respond  to  any  tempo¬ 
ral  variations  in  the  channel,  this  open-loop  control  algorithm  probably  would  not 
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perform  well  in  many  ocean  environments. 

One  of  the  best  experiments  for  single-mode  generation  in  a  laboratory  setting 
was  performed  by  Gazanhes  and  Gamier,  [16].  They  used  a  tank  with  only  57mm  of 
water  over  a  sand  bottom.  At  their  working  frequency  of  124  kHz,  this  environment 
only  supported  5  trapped  modes,  making  control  and  estimation  of  the  modes  much 
more  tractable  than  the  Glay  and  Huang  experiment.  Using  the  Pekeris  model  of 
propagation,  Gazanhes  and  Gamier  computed  the  mode  shapes  for  the  channel  and 
used  samples  of  these  shapes  to  weight  15  independent  piezoelectric  transducers. 
Although  this  is  not  the  optimal  weighting  in  the  least-squares  sense,  it  is  a  much 
better  discrete  approximation  to  the  orthogonality  condition  than  Clay  and  Huang’s 
very  coarse  sampling  and  weighting.  Gazanhes  and  Gamier  confirmed  that  they 
excited  only  a  single  mode  using  synthetic  aperture  arrays  in  both  range  and  depth. 
Although  their  control  scheme  is  also  open  loop,  it  works  well  because  they  have  a 
known  time-invariant  environment.  This  is  encouraging  evidence  that  it  is  possible 
to  excite  a  single  mode  given  an  accurate  channel  estimate. 

This  concludes  our  review  of  previous  work  on  mode  propagation  both  in  the 
ocean  and  in  laboratory  settings.  The  next  section  describes  the  Green’s  function 
for  underwater  acoustics.  The  Green’s  function  is  the  foundation  of  our  model  for 
the  channel  between  the  source  array  and  feedback  array. 

1.2  Green’s  Function 

The  Green’s  function  characterizes  the  acoustic  propagation  between  a  source  and 
receiver  at  a  given  frequency.  To  find  this  function  in  regions  of  constant  density, 
simplify  the  Helmholtz  equation  (Eq.  1.1)  to  obtain 

-k  A:^(r)]  p(r)  = -47r5(r  -  To)  (1-12) 
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for  a  point  source  at  ro-  The  solution  to  this  equation,  p(r)  =  (?(r,ro),  is  the 
Green’s  function  [22],  [32].  In  linear  system  theory  terms,  the  Green’s  function  is 
the  transfer  function  between  a  point  source  at  Tq  and  a  receiver  at  r  evaluated  at 
a  single  frequency  w.  The  Green’s  function  is  a  very  general  model  for  propagation, 
and  is  still  valid  in  many  environments  where  a  discrete  set  of  normal  modes  cannot 
accurately  model  propagation.  For  a  driving  term  / (ro)  it  can  be  shown  the  pressure 
field  is 


^  Iv  ^o)f(ro)dro  +  ^  /  fecr,  ro)^|^  -  p(ro)^^^ 

Jvo  47r  Jso[  dno  dno 


dro.  (1.13) 


The  first  integral  incorporates  the  contribution  of  all  sources  in  the  volume  in  ques¬ 
tion  Vo,  while  the  second  integral  includes  the  effect  of  the  boundary  conditions  on 
-So,  the  surface  surrounding  Vq. 

If  the  acoustic  channel  between  source  and  receiver  arrays  is  a  linear  system  with 
respect  to  the  complex  weight  vector  of  the  source  array,  the  contribution  of  the 
second  integral  in  Eq.  1.13  must  be  zero,  and  there  must  not  be  any  other  significant 
sources  in  the  volume.  The  condition  that  the  second  integral  of  Eq.  1.13  be  zero  is 
the  spatial  analog  of  the  initial  rest  boundary  condition  for  a  time-domain  linear  sys¬ 
tem.  In  physical  terms,  this  corresponds  to  an  absence  of  sources  outside  the  volume 
in  question,  Vq.  In  general,  the  ocean  may  be  considered  time-invariant  for  the  rela¬ 
tively  brief  propagation  times  between  the  feedback  and  source  arrays.  Under  these 
conditions,  the  channel  is  well-modeled  by  a  linear,  time-invariant  (LTI)  system. 
Complex  exponentials  are  eigenfunctions  of  LTI  systems,  so  the  Green’s  function 
completely  characterizes  the  behavior  of  the  channel  at  the  frequency  of  excitation. 
Because  it  is  a  very  general  model  for  propagation,  the  Green’s  function  can  incorpo¬ 
rate  effects  from  any  spatial  frequencies  in  the  modal  continuum  which  are  not  one  of 
the  trapped  modes,  but  still  which  couple  energy  back  into  the  propagating  modes, 
in  addition  to  summarizing  any  coupling  of  energy  among  the  propagating  modes. 
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In  environments  which  are  completely  represented  by  a  discrete  set  of  normal  modes, 
the  depth-dependent  Green’s  function  has  poles  at  the  horizontal  wavenumbers  of 
the  modes,  indicating  the  pressure  field  consists  predominantly  of  energy  propagat¬ 
ing  at  those  spatial  frequencies.  Because  the  set  of  environments  where  the  Green’s 
function  accurately  models  propagation  is  a  superset  of  those  accurately  modeled 
by  either  the  adiabatic  or  coupled  mode  models,  it  is  a  good  model  for  the  control 
algorithm  to  use  to  summarize  the  acoustic  propagation  through  the  feedback  vol¬ 
ume.  The  model  allows  single  mode  excitation  in  any  environment  well-modeled  by 
either  adiabatic  or  coupled  mode  propagation,  and  possibly  in  some  more  complex 
environments  as  well. 


1.3  Mode  Filters 


The  problem  of  estimating  the  coefficients  of  the  modes  propagating  at  a  given  lo¬ 
cation  from  samples  of  the  pressure  held  obtained  at  hydrophones,  known  as  mode 
filtering,  is  a  common  one  in  ocean  acoustics.  The  feedback  control  algorithm  must 
determine  how  the  error  between  the  desired  and  observed  pressure  profile  at  the 
feedback  array  is  related  to  the  error  between  the  desired  and  excited  mode  coeffi¬ 
cients.  These  two  problems  are  closely  related,  and  in  addition  to  reviewing  common 
algorithms  used  for  mode  filtering,  this  section  will  also  examine  the  similarity  be¬ 
tween  these  problems. 

The  pressure  field  generated  by  M  propagating  modes  spatially  sampled  by  a 
vertical  array  of  N  hydrophones  can  be  written  as 


p(^l) 

di 

n{zi) 

: 

: 

+ 

: 

pM 

d-M 

n(zAr) 

(1.14) 
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or  in  vector  notation 


p  =  +  n,  (1.15) 

where  n  is  the  vector  of  observation  noise  at  the  hydrophone  locations. 

Linear  mode  filters  estimate  the  mode  coefficients  as  a  linear  function  of  the 
observed  pressure  samples.  For  Sections  1.3.1  through  1.3.3,  the  mode  coefficient 
vector  will  be  considered  a  deterministic  but  unknown  quantity  to  be  estimated. 
In  This  linear  function  can  be  represented  by  a  matrix  H  multiplying  the  observed 
pressure  field  p,  i.e., 

d  =  Hp  =  H^d  +  Hn.  (1.16) 

The  various  mode  filters  discussed  in  this  section  correspond  to  different  choices  for  H 
proposed  in  the  literature  of  generalized  inverses  [33],  [34]  and  stochastic  estimation 
[35].  For  mode  filters  of  this  form,  the  covariance  of  the  mode  coefficient  vector 
estimate  =  F'ldd^}  -  F'{d}£'{d”}  is 


Kaa  =  (1.17) 

where  is  the  spatial  covariance  of  the  noise  vector  n,  and  (•)^  denotes  the 
Hermitian,  or  conjugate-transpose,  operator.  Thus,  the  covariance  of  the  estimated 
mode  coefficients  depends  entirely  on  the  noise  process  covariance  Kn„  and  the  mode 
filter  H.  The  remainder  of  this  section  reviews  the  common  choices  for  H  in  the  ocean 
acoustics  literature. 

1.3.1  Sampled  Mode  Shapes  Mode  Filter 

A  common  choice  for  H  in  mode  propagation  experiments  is  The  motivation  for 
this  choice  is  that  the  mode  shapes  are  orthonormal  when  integrated  as  a  function 
of  continuous  depth  and  weighted  by  the  density  profile  [22].  However,  spatially 
sampling  the  modes  does  not  necessarily  preserve  the  orthogonality  of  the  modes. 
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i.e.,  ^  I.  This  lack  of  orthogonality  can  introduce  a  bias  into  the  mode 

estimate.  Assuming  the  noise  vector  n  is  zero-mean,  the  expected  value  of  the  mode 
coefficient  estimator  is 

E{d}  =  (1.18) 

giving  a  bias  of  (/—  Theoretically,  as  more  hydrophones  are  added  to  sam¬ 

ple  the  water  column  more  finely  in  depth,  the  sampled  modes  matrix  ^  becomes 
arbitrarily  close  to  orthogonal,  making  the  effect  of  the  bias  negligible.  Realistically, 
it  is  impractical  to  deploy  a  vertical  array  of  hydrophones  spanning  the  entire  sed¬ 
iment  layer.  Thus,  even  in  the  limiting  case  of  a  continuous  array  of  hydrophones 
spanning  the  water  column,  some  bias  will  still  be  present  due  to  the  unsampled 
pressure  field  in  the  bottom  [36].  Many  experiments  including  Ferris  [8],  and  Clay 
and  Huang  [15]  used  the  sampled  mode  shape  mode  filter,  assuming  the  samples 
of  the  orthogonal  mode  functions  are  themselves  orthogonal  without  examining  the 
potential  bias  of  the  sampled  mode  shape  filter.  Tindle,  Hobaek,  and  Muir  [31]  and 
Gazanhes  and  Gamier  [16]  used  the  same  mode  filter  but  examined  the  “cross-talk” 
introduced  by  sampling  to  verify  the  bias  was  negligible  for  the  purposes  of  their 
experiments. 

The  mode  coefficients  may  be  difficult  to  determine  from  the  pressure  samples 
taken  at  the  hydrophone  locations  if  the  hydrophone  array  either  has  fewer  hy¬ 
drophones  than  the  number  of  propagating  modes,  or  if  the  hydrophones  do  not 
fully  span  the  water  depth.  The  first  difficulty  will  be  referred  to  as  undersampling, 
while  the  second  will  be  referred  to  as  poorly-conditioned  sampling.  The  singular 
value  decomposition  (SVD)  [37],  an  orthogonal  matrix  factorization,  provides  in¬ 
sight  on  how  these  problems  affect  the  mode  estimator.  The  SVD  factors  ^  as 
where  and  are  orthonormal  matrices.  Assuming  the  number 


28 


of  hydrophones  exceeds  the  number  of  trapped  modes,  i.e.,  N  >  M, 


(1.19) 


where  cr^i  >  <7^2  >...>  >  0  are  called  the  singular  values  of  In  the 

undersampled  case,  with  N  <  M,  the  mode  filtering  problem  is  underdetermined  and 
d  cannot  be  uniquely  determined  from  p.  Practically,  this  is  usually  not  a  problem 
for  mid-to-low  frequencies  in  the  shallow  water  environment,  since  it  is  not  difficult 
or  expensive  to  construct  an  array  with  more  hydrophones  than  the  environment 
has  propagating  modes.  Even  if  >  M,  the  choice  of  hydrophone  locations  may 
give  a  poorly-conditioned  problem  for  determining  d.  Although  the  problem  is  now 
nominally  over-determined,  poorly-conditioned  sampling  will  cause  to  become 
rank-deficient  or  nearly  rank-deficient.  This  will  make  it  impossible  to  estimate  d 
reliably  in  the  presence  of  noise.  Golub  and  Van  Loan  [37]  demonstrate  poorly- 
conditioned  sampling  results  in  one  or  more  of  the  singular  values  am  becoming  very 
small  or  zero.  Rewriting  the  bias  R(d)  as 
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The  covariance  of  the  sampled  mode  shapes  estimator  depends  on  the  covariance 
of  the  noise  vector  Knn-  Two  common  models  for  noise  in  the  shallow  water  channel 
are  spatially  white  noise  {K^n  =  <7^I)  and  the  Kuperman-Ingenito  surface  noise 
model  [6].  The  Kuperman-Ingenito  model  proposes  that  when  the  noise  is  dominated 
by  surface  generated  noise,  i.e.,  noise  due  to  wind,  surface  waves,  etc. . . ,  it  will  have 
a  spatial  covariance  of  the  form 


(t|  0 

d\ 


0  (T? 

0,2 


0  cr? 

OM 


(1.21) 


where  cr|^ , . . . ,  are  functions  of  the  mode  profiles  and  surface  noise  processes. 
The  vector  d  is  the  mode  coefficients  of  the  noise  process  at  the  array. 

For  the  spatially  white  case,  Eq.  1.17  becomes 


M 

^dd  =  (1-22) 

m=l 

where  are  the  columns  of  V^.  The  Kuperman-Ingenito  model  gives 


^dd  = 


M  \ 

71=1  / 


0  ...  0  (T? 


M  \ 

)  •  (1-23) 

71=1  / 


Poorly-conditioned  sampling  decreases  some  values  of  cr^m,  but  its  effect  on  the  es¬ 
timator  covariance  does  not  grow  more  extreme  as  the  conditioning  of  'i'  becomes 


more  severe.  Once  a  particular  singular  value  grows  negligible  compared  to  other 
singular  values,  further  decreases  in  its  value  are  insignificant  in  the  covariance  com¬ 
putations  in  Eqs.  1.22  and  1.23.  At  the  other  extreme,  when  the  channel  is  highly 
oversampled,  i.e.,  N  oo,  the  approach  1,  so 

M 

^  =  I  (1.24) 

m=l 

and  Eqs.  1.22  and  1.23  simplify  to  and  respectively. 

Thus,  poorly-conditioned  sampling  of  the  water  column  by  the  hydrophone  array 
will  increase  the  bias  of  the  sampled  mode  shape  filter,  but  will  not  effect  the  estima¬ 
tor  covariance  significantly.  Oversampling  the  channel  so  that  the  inner  product  of 
the  sampled  mode  shapes  closely  approximates  the  integral  of  the  continuous  mode 
shapes  can  reduce  the  bias  to  leave  only  a  small  contribution  from  the  unsampled 
bottom,  and  result  in  an  estimator  covariance  reflecting  the  underlying  nature  of  the 
noise  process. 

1.3.2  Pseudo-Inverse  Mode  Filter 

The  pseudo-inverse  mode  filter  results  from  choosing  d  to  minimize  the  squared 
error  between  "I'd  and  p.  Tindle  et  al.  [9]  appears  to  be  the  first  reference  in  the 
ocean  acoustics  literature  to  formulate  the  mode  estimation  problem  in  this  least 
squares  sense.  The  resulting  mode  filter  H  =  denoted  is  called 

the  pseudo-inverse  or  Penrose-Moore  inverse  of  ^  [37],  [33].  This  name  results  from 
the  fact  =  I.  A  consequence  of  this  property  is  that  if  E  {n}  =  0,  the  mode 
estimate  obtained  using  H  =  ’4^1  is  unbiased.  However,  if  the  hydrophone  array  gives 
a  poorly-conditioned  sampling  of  the  modes,  this  estimator  becomes  very  sensitive 
to  noise. 
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For  the  spatially  white  noise  model,  Eq.  1.17  gives 


^dd  = 

M 

=  (1.25) 

m=l 

Although  this  estimate  will  remain  unbiased  so  long  as  N  >  M,  for  very  short 
aperture  arrays  the  estimator  will  have  a  very  large  covariance  as  one  more  or  more 
singular  values  aq,„i  go  to  zero.  When  a  full  aperture  array  oversamples  the  channel, 
the  singular  values  of  approach  1,  and  Eq.  1.25  approaches  cr^I. 

Substituting  the  Kuperman-Ingenito  noise  model  and  the  pseudo-inverse  mode 
filter  into  Eq.  1.17  gives 


=  Kii,  (1.26) 

which  is  intuitively  sensible,  as  is  the  covariance  of  the  noise  process  as  it  is 
coupled  into  the  channel  by  the  modes.  Ideally,  this  covariance  is  unchanged  by 
reductions  in  the  array  aperture.  Practically,  the  mode  filter  is  usually  based  on 
an  estimate  of  computed  by  numerical  integration  of  an  observed  or  estimated 
sound  speed  profile.  As  the  <t„,’s  of  the  true  ^  grow  smaller  as  the  array  aperture 
decreases  and  the  sampling  becomes  poorly-conditioned,  the  bias  and  covariance  of 
the  pseudo-inverse  mode  filter  can  become  very  sensitive  to  errors  between  the 
obtained  by  numerical  integration  and  the  actual  ^  of  the  ocean  channel. 

For  the  spatially  white  Gaussian  noise  case,  the  pseudo-inverse  mode  filter  can  be 
shown  to  be  the  maximum-likelihood  (ML)  estimator.  In  general,  it  is  not  an  efficient 
estimator  and  does  not  attain  the  Cramer-Rao  Bound  (CRB)  for  the  variance  of 
unbiased  estimators  unless  all  the  singular  values  are  equal.  Poorly-conditioned 
sampling  of  the  channel  can  push  the  covariance  of  this  mode  filter  arbitrarily  far 
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from  the  CRB  as  some  of  the  singular  values  grow  very  small.  If  an  efficient  estimator 
exists,  it  is  the  ML  estimator  [35];  since  the  ML  estimator  is  not  efficient,  no  efficient 
estimator  exists  for  the  spatially  white  noise  model. 

The  pseudo-inverse  mode  filter  is  also  the  ML  estimator  for  the  Kuperman- 
Ingenito  noise  model.  For  this  case,  this  conclusion  is  not  very  informative  since 
n  is  defined  to  be  in  the  range  of  ^  from  the  formulation  of  the  noise  model.  The 
existence  of  the  ML  estimator  depends  on  the  Kuperman-Ingenito  model  perfectly 
describing  the  noise  process,  since  if  n  contains  any  component  in  the  orthogonal 
complement  to  the  range  of  \E^,  the  conditional  probability  density  pp|D(p|d)  =  0 
[38],  and  the  ML  estimate  is  meaningless  since  no  set  of  mode  coefficients  d  could 
have  produced  the  observed  signal.  For  this  noise  model,  the  Fisher  Information 
Matrix  can  be  shown  to  be  so  the  CRB  for  the  mode  estimator  is  the 
variance  of  the  underlying  noise  process  in  each  mode.  The  pseudo-inverse  mode 
filter  attains  this  bound  on  the  variance,  and  is  an  efficient  estimator. 

1.3.3  Diagonal  Weighting 

As  noted  above,  when  the  array  has  poorly-conditioned  sampling  of  the  modes,  one 
or  more  of  the  singular  values  of  ^  will  grow  very  small.  This  causes  to  be 

singular  or  nearly  singular,  and  the  computation  of  the  inverse  of  this  matrix  can  grow 
numerically  sensitive.  One  method  of  compensating  for  this  is  to  modify  the  error 
function  being  minimized  to  include  a  term  proportional  to  the  magnitude  squared 
of  the  estimated  mode  coefficient  vector  d  [39].  The  quantity  to  be  minimized  is 
then 

e  =  ]|p-^dlp-h/?]|dl|2  (1.27) 
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where  /?  is  a  scale  factor  indicative  of  the  relative  importance  of  the  two  terms  in 
the  error  expression.  The  estimator  minimizing  this  quantity  is 

dow  =  +  I3iy^  (1.28) 

This  expression  is  very  similar  to  the  pseudo-inverse,  except  for  a  small  diagonal  ma¬ 
trix  pi  has  been  added  to  before  inversion  to  alleviate  conditioning  problems. 
The  pi  term  is  often  referred  to  as  the  white  noise  sensitivity  term.  The  addition 
of  this  term  places  a  lower  bound  of  p  on  the  singular  values  of  -f  pi),  and 

thus  limits  the  covariance  of  the  estimator  shown  in  Eq.  1.25,  since  no  for  the 
diagonally-weighted  inverse  can  exceed  P^"^.  For  this  reason,  this  approach  is  often 
referred  to  as  diagonal  loading  or  weighting.  While  this  estimator  does  not  possess 
many  of  the  nice  theoretical  properties  of  the  pseudo-inverse  mode  filter,  in  many 
scenarios  it  is  computationally  more  stable,  especially  at  relatively  high  noise  levels. 

1.3.4  Maximum  A  Posteriori  Mode  Filters 

The  maximum  a  posteriori  (MAP)  mode  filter  chooses  dMAP  to  maximize  the  prob¬ 
ability  density  function  pD|p(d|p)  based  on  the  observed  pressure  field,  p.  This 
approach  assumes  the  mode  coefficient  vector  d  is  an  unknown,  random  quantity  to 
be  determined.  The  mode  filter  also  assumes  knowledge  of  the  probability  density 
functions  of  both  the  stochastic  processes  generating  the  mode  coefficients  d  and  the 
noise  n.  For  the  case  when  both  these  processes  are  Gaussian,  the  MAP  estimate 
is  identical  to  the  minimum  mean-squared  error  estimate  (MMSE)  [35].  A  very  im¬ 
portant  feature  of  the  MAP  mode  filter  is  that  in  many  environments  it  gracefully 
transitions  from  the  pseudo-inverse  filter  to  the  sampled  mode  shape  filter  as  the  con¬ 
ditioning  of  4^  deteriorates.  Between  these  extremes,  the  performance  of  the  MAP 
filter  exceeds  either  of  the  other  two  filters,  and  the  MAP  filter  never  exhibits  the 
poor  performance  shown  by  the  pseudo-inverse  or  sampled  mode  shape  filters  when 
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these  filters  are  applied  to  a  problem  whose  conditioning  is  inappropriate  for  their 
strengths.  If  the  mode  coefficients  are  well-modeled  by  a  Gaussian  process  with  zero 
mean  and  covariance  K^di  while  the  noise  is  also  Gaussian,  zero-mean,  uncorrelated 
with  d,  and  has  covariance  the  MAP  mode  filter  can  be  written  as 

^MAP  =  (1.29) 

where 

(1-30) 

Some  insight  into  the  performance  of  this  mode  filter  may  be  gained  by  consider¬ 
ing  the  somewhat  unrealistic  case  when  the  modes  are  independent  and  identically 
distributed,  i.e.,  K^d  =  o’jl,  and  the  noise  is  spatially  white  with  K^n  =  cr^I-  As¬ 
suming  there  are  more  hydrophones  than  modes  {M  >  N),  Eq.  1.29  reduces  to 


dMAP  = 


Eq.  1.31  has  an  appealing  interpretation  as  a  generalization  of  the  discrete  spatial 
Wiener  filter  (DSWF)  [40].  Multiplying  by  rotates  the  problem  into  the  coordi¬ 
nate  frame  where  the  spatial  components  are  uncorrelated.  Each  component  is  then 
weighted  by  the  Wiener  gain  for  the  ratio  of  the  mode  power  to  the  noise  power  for 
that  component  <^%rn^d/ +  <^n)-  These  estimates  of  the  components  are  then 
multiplied  by  the  inverse  singular  values  before  being  transformed  into  mode 
coefficients  by  V^.  For  the  case  when  all  the  singular  values  are  1,  =  I  and 

is  the  appropriate  set  of  samples  of  complex  exponentials,  Eq.  1.31  reduces  exactly 
to  the  DSWF. 


-2^2 


0 
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U^“p.  (1.31) 
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Eq.  1.31  provides  theoretical  justification  for  the  pragmatically-motivated  mode 
filtering  algorithm  proposed  by  T.  C.  Yang  in  [41].  Yang’s  algorithm  suggested 
deleting  very  small  eigenvalues  of  before  inverting  this  product  for  the  pseudo¬ 
inverse  mode  filter.  His  paper  proposes  this  approach  to  ensure  numerical  stability  of 
the  inverse  and  uses  some  rough  a  priori  assumptions  about  the  values  of  the  mode 
coefficients.  His  approach  is  a  limiting  case  of  Eq.  1.31  as  some  of  the  singular  values 
of  get  very  small  compared  to  and  cr„.  Under  these  conditions,  Yang’s  ad  hoc 
“dropped  eigenvalue”  method  is  an  approximation  to  the  MAP  mode  filter  obtained 
when  the  array  gives  a  poorly-conditioned  sampling  of  the  mode  shapes  and  some 
of  the  diagonal  elements  in  Eq.  1.31  approach  zero. 

If  the  Kuperman-Ingenito  noise  model  applies,  and  the  process  generating  the 
modes  has  diagonal  covariance  with  elements  , . . . ,  ,  the  MAP  estimator 

is 


dMAP  = 


.+0-1 


+crj 
<^2  d2 


J 


(1.32) 


Because  both  the  signal  and  noise  processes  are  independent  random  variables  trans¬ 
formed  by  the  same  linear  transformation  (^),  multiplying  p  by  the  pseudo-inverse 
performs  the  discrete  spatial  Karhunen-Loeve  transform,  decorrelating  the  mode  es¬ 
timation  problem  into  M  independent  problems.  Unlike  the  spatially  white  noise 
model,  the  physical  basis  of  interest  (mode  space)  coincides  with  the  mathematical 
basis  in  which  the  underlying  processes  are  independent.  The  diagonal  elements  in 
Eq.  1.32  become  the  Wiener  gains  for  each  of  the  independent  components.  Thus 
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the  structure  of  Eq.  1.32  shows  that  the  estimation  problem  becomes  M  indepen¬ 
dent  MAP  problems  once  p  is  multiplied  by  When  the  random  processes  have 
probability  density  functions  which  are  symmetric  about  their  means,  such  as  the 
Gaussian  distribution,  Eqs.  1.31  and  1.32  are  also  the  minimum  mean-squared  error 
(MMSE)  solutions  [35]. 

1.3.5  Relation  of  Mode  Filters  to  Target  Pressure  Vector 

As  mentioned  earlier,  the  variety  of  mode  filtering  algorithms  corresponds  to  a  variety 
of  choices  for  the  target  pressure  vector  for  the  control  algorithm.  It  is  possible  to 
formulate  the  control  problem  in  terms  of  the  desired  pressure  field  rather  than 
the  desired  mode  coefficients.  Section  2.2  discusses  the  circumstances  when  it  is 
preferable  to  use  this  approach.  The  control  algorithm  uses  an  estimate  of  how  the 
pressure  field  observed  at  the  hydrophone  p  is  a  function  of  the  source  array  weight 
vector  w  to  choose  w  to  obtain  some  target  pressure  vector  p^.  Formulating  the 
problem  as  a  mode  filtering  problem  at  the  feedback  array  gives  insight  into  the  best 
choice  of  p</.  The  target  pressure  vector  should  be  the  pressure  vector  which  would 
yield  the  desired  mode  coefficient  vector  when  the  most  appropriate  mode  filter 
for  the  ocean  environment  and  feedback  array  geometry  is  applied. 

The  most  obvious  choice  for  p^  is  This  target  pressure  vector  corre¬ 
sponds  to  the  pseudo-inverse  mode  filter,  since  =  ^f^d^  =  d^.  Thus,  if 

this  choice  for  pd  were  observed  at  the  feedback  array,  the  pseudo-inverse  mode  filter 
would  indicate  the  pressure  consisted  of  the  desired  modes.  Alternatively,  the  choice 
Pd  =  corresponds  to  the  sampled  mode  shape  filter,  since  would 

give  =  dd-  Similar  expressions  can  be  derived  for  the  MAP  mode 

filters  with  the  various  noise  models. 

The  algorithms  described  in  this  thesis  will  use  Pd  =  ^dd-  The  feedback  array 
is  assumed  to  sample  the  channel  sufficiently  that  the  pseudo-inverse  mode  filter 
is  the  best  choice  of  mode  filter.  This  is  a  reasonable  assumption  for  the  relatively 
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shallow  water  environments  and  mid-frequency  transmissions  examined  in  this  thesis. 
Also,  over  the  relatively  short  range  between  the  source  array  and  feedback  array, 
the  signal-to-noise  (SNR)  ratio  at  the  feedback  array  should  be  high  if  the  signal 
is  going  to  be  observable  at  the  distant  observation  array  as  shown  in  Fig.  1-1.  At 
high  SNR,  the  MAP  formulations  reduce  asymptotically  to  If  these 

favorable  conditions  do  not  apply,  the  target  pressure  vector  should  be  chosen  based 
upon  the  appropriate  mode  filtering  algorithm  for  the  environment  and  geometry, 
but  generally  this  means  the  algorithm  will  not  be  able  to  produce  a  single  mode 
with  high  fidelity. 


1.4  Control  Theory 

The  field  of  control  theory  provides  some  important  insights  into  the  single  mode 
excitation  problem.  The  classification  of  control  systems  as  either  open-loop  and 
closed-loop  systems  is  one  of  the  fundamental  distinctions  in  control  theory  [42]. 
Both  kinds  of  systems  use  a  control  law  for  determining  the  control  input  to  attempt 
to  attain  the  desired  behavior  of  the  output  of  the  plant,  or  system  controlled.  For 
the  mode  excitation  problem,  the  input  is  the  source  array  weights,  w,  the  plant  is 
the  ocean  between  the  source  array  and  the  feedback  array,  the  output  is  the  samples 
of  the  pressure  field  observed  at  the  feedback  array,  p,  and  the  desired  behavior  is 
the  target  pressure  p<i.  For  an  open-loop  system,  the  input  is  chosen  based  solely 
on  a  priori  or  assumed  knowledge  of  the  system.  All  prior  work  on  single  mode 
excitation  used  open-loop  control  [15],  [16].  The  behavior  of  a  plant  controlled  with 
this  approach  can  be  very  sensitive  to  inaccuracies  in  this  a  priori  information  or  in 
the  model  of  the  plant.  The  algorithm  proposed  in  this  thesis  uses  feedback  control 
to  obtain  the  desired  performance.  Feedback  control  modifies  the  input  to  the  plant 
based  on  the  difference  between  the  observed  and  desired  output  of  the  plant.  This 
allows  the  control  system  to  compensate  for  inaccuracies  in  the  a  priori  information 
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about  the  plant  as  well  as  respond  to  time-variations  in  the  plant.  Figs.  1-2  and  1-3 
below  illustrate  the  two  approaches  for  the  single  mode  excitation  problem. 


Figure  1-2;  Open-loop  Control  For  Single  Mode  Excitation 


Figure  1-3;  Closed-loop  Control  For  Single  Mode  Excitation 

Another  useful  concept  from  control  theory  is  the  use  of  state  space  models  for 
systems  with  complicated  but  linear  dynamics.  Many  feedback  control  algorithms 
use  such  a  model  to  estimate  the  behavior  of  a  plant  they  wish  to  control.  Applying 
this  model  to  the  acoustic  propagation  between  the  source  and  feedback  arrays  pro¬ 
vides  important  insight  into  several  aspects  of  the  single  mode  excitation  problem. 
For  the  discrete-time  case,  the  state-space  model  describes  the  evolution  of  the  1  x  M 
state  vector  x[n]  by 

x[n]  =  A[n  -  l]x[n  -  1]  +  B[n]u[n]  -I-  f[n|  (1.33) 

where  u[n]  is  the  1  x  L  input  vector  for  the  system,  and  f[n]  is  a  zero-mean  white 
Gaussian  noise  with  covariance  Pff  which  is  uncorrelated  with  all  other  signals  in 
the  problem.  The  output  of  the  system  y[n]  is  a  1  x  AT  vector 

y[n]  =  C[n]x[n]  +  D[n]u[n]  +  n[n],  (1.34) 
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where  n[n]  is  zero-mean,  white  Gaussian  observation  noise  uncorrelated  with  all 
other  signals  in  the  problem,  with  covariance  P„„.  Defining  the  state  vector  x[n]  as 
the  modes  d[n]  present  at  the  feedback  array  at  the  iteration,  the  input  as  w[n], 
the  complex  source  weights,  and  the  output  y[n]  as  p[n],  the  pressure  observed  at 
the  feedback  array,  the  state  space  model  for  the  modal  propagation  between  the 
source  and  feedback  arrays  simplifies  to 

d[n]  =  B[n]w[n]  +  f[n]  (1.35) 

p[n]  =  ^[n]d[n]  +  n[n].  (1.36) 

In  these  equations,  the  B[n]w[n]  term  is  the  portion  of  the  total  pressure  field  excited 
the  source  array,  the  f[n]  term  is  the  propagating  background  noise  which  satisfies 
the  Kuperman-Ingenito  noise  model,  and  n[n]  represents  the  sensor  noise.  The  devel¬ 
opment  that  follows  does  not  fully  exploit  this  general  model,  but  it  is  a  conceptually 
appealing  model  for  scenarios  when  there  are  different  noise  sources  with  differing 
spatial  covariance  structures.  The  matrix  B[n]  incorporates  the  effects  of  the  mode 
shapes  at  the  source  array  on  the  energy  initially  excited  in  each  mode  (Eq.  1.7), 
as  well  as  any  coupling  of  energy  among  the  discrete  modes  between  the  source  and 
feedback  arrays.  Note  that  Eq.  1.36  has  the  same  form  as  the  mode  filtering  problem 
(Eq.  1.15).  Modeling  the  modal  propagation  in  this  way  allows  the  application  of 
several  powerful  results  from  control  theory  to  the  single  mode  excitation  problem. 

Reachability  and  observability  are  two  control  theory  concepts  germane  to  the 
mode  excitation  problem  as  formulated  above  in  Eqs.  1.35  and  1.36.  A  system  is 
said  to  be  reachable  over  the  interval  [no,  nj  if  it  is  possible  to  drive  the  system  to 
any  final  state  x[ni]  from  any  initial  state  x[no]  using  u[n].  The  effect  of  the  input 
u[n]  on  the  state  x[ni]  can  be  written  as  ^[ni,n]B[n]u[n],  where  #[ni,n]  is  the 
state-transition  matrix  from  time  n  to  rii  for  the  undriven  system.  In  the  absence 
of  inputs,  x[ni]  =  $[ni,  n]x[n].  The  total  effect  of  the  input  over  the  interval  [no,  ni] 


40 


on  the  state  at  ni  can  be  written  as 


x[^i]  —  B[ni]  <&[ni,ni  —  l]B[ni  -  1]  ...  #[ni,no]B[no] 

(1.37) 

For  a  system  to  be  reachable,  the  partitioned  matrix  containing  the  B  and  $  terms 
must  have  full  rank.  The  controllability  gramian  Wi[no,  ni\  of  a  system  is  defined  to 
be  the  correlation  of  this  matrix.  A  system  is  reachable  if  and  only  if  Wi[no,  rii]  has 
full  rank.  For  the  simple  state  space  model  of  Eqs.  1.35  and  1.36,  the  controllability 
gramian  is 

Wi[no,ni]  =  BfnijB'^fni],  (1.38) 

so  the  system  is  reachable  over  [no,ni]  if  and  only  if  B[ni]  has  full  row  rank  M.  A 
system  is  defined  to  be  completely  reachable  if  it  is  reachable  over  some  interval, 
and  is  uniformly  completely  reachable  if  there  exists  some  a  and  m  such  that  the 
smallest  eigenvalue  of  Wi[n,  n  +  m]  is  greater  than  a  for  all  n.  Because  the  modes 
present  at  the  feedback  array  for  this  model  depend  only  on  the  input  weights  from 
the  current  iteration,  properties  like  uniform  complete  reachability  depend  solely  on 

BH- 

A  system  is  described  as  observable  over  the  interval  [no,ni]  if  the  state  x[no] 
can  be  uniquely  determined  from  knowing  u[n]  and  y[n]  over  this  interval.  It  can 
be  shown  this  property  is  equivalent  to  the  observability  gramian  M[no,ni]  having 
full  rank  M  [43].  The  very  simple  state  space  model  of  mode  propagation  has  the 
following  observability  gramian: 

M[no,ni]  =  ^“[no]^[no]. 
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Thus,  the  coefficients  of  the  modes  present  at  the  feedback  array  can  only  be  uniquely 
determined  if  the  sampled  mode  shape  matrix  has  full  column  rank. 

Framing  the  mode  excitation  problem  in  the  state  space  model  clearly  separates 
the  two  primary  issues  of  the  problem.  The  first  issue  is  whether  it  is  possible  to  excite 
the  desired  modes  in  the  channel.  The  second  issue  is  whether  the  control  algorithm 
can  recognize  the  desired  mode  profile  by  observing  samples  of  the  pressure  field. 
The  first  issue  corresponds  to  the  reachability  of  the  system  as  described  by  the  state 
space  model  for  mode  propagation.  The  reachability  depends  on  the  source  array 
geometry  and  the  propagation  conditions  between  the  source  array  and  feedback 
array.  For  the  ranges  and  depths  discussed  in  this  thesis,  the  propagation  is  usually 
predominantly  determined  by  the  bottom  composition  and  bathymetry.  Since  the 
source  array  is  the  only  factor  effecting  reachability  under  the  control  of  the  scientist, 
care  must  be  taken  to  design  the  source  array  to  ensure  reachability  over  as  wide  a 
range  of  ocean  conditions  as  possible. 

The  ability  of  the  algorithm  to  recognize  the  desired  mode  profile  from  pressure 
observations  corresponds  to  the  observability  of  the  system  in  the  state  space  formu¬ 
lation.  The  observability  of  the  system  depends  on  the  mode  shapes  and  feedback 
array  geometry.  The  mode  shapes  are  a  function  of  the  sound  speed  profile  and  bot¬ 
tom  properties  at  the  array  location,  and  are  not  under  the  control  of  the  observer. 
Thus,  the  feedback  array  must  be  designed  carefully  to  guarantee  observability  over 
as  wide  a  range  of  sound  speed  profiles  and  bottom  properties  as  possible  to  ensure 
the  experiment’s  success. 

An  informative  experiment  and  persistency  of  input  are  two  additional  important 
concepts  from  system  identification  and  control  theory.  An  experiment  is  said  to  be 
informative  if  it  allows  the  unique  identification  of  a  single  set  of  system  model 
parameters  out  of  the  class  of  proposed  system  models  [44].  A  necessary  condition 
for  an  experiment  to  be  informative  is  that  the  control  input  sequence  be  persistent. 
A  persistent  input  sequence  is  one  whose  autocorrelation  matrix  has  full  rank.  This 
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criteria  will  be  formulated  more  precisely  in  Chapter  2  where  the  class  of  channel 
models  used  by  the  control  algorithm  is  presented.  The  conflicting  requirements  of 
keeping  the  source  weights  fixed  to  obtain  a  high  fidelity  mode  and  varying  them 
enough  to  get  a  persistent  input  for  system  identification  will  be  an  important  issue 
for  the  control  algorithm  to  address. 

1.5  Matched  Signals 

Single  mode  excitation  at  a  feedback  array  location  can  be  viewed  as  a  general¬ 
ization  from  the  time  domain  to  the  mode  domain  of  the  matched  signal  problem 
addressed  by  Parvulescu  and  Clay  [45],  [46],  [47].  Parvulescu  and  Clay  investigated 
the  time-domain  dispersion  due  to  propagation  through  underwater  acoustic  chan¬ 
nels.  Specifically,  they  were  interested  in  finding  a  signal  which  when  transmitted 
from  a  single  source  through  a  given  channel  would  give  an  “impulse-like”  signal  at 
a  single  receiving  hydrophone.  They  observed  that  when  they  transmitted  a  short 
impulsive  waveform,  the  received  signal  exhibited  the  well-known  multi-path  arrival 
structure.  The  matched  filter  whose  impulse  response  is  the  time-reversed  version 
of  the  received  signal  is  the  best  linear  detector  for  this  signal  in  the  presence  of 
white  noise  [35].  Parvulescu  and  Clay  noted  that  if  the  transmitted  signal  is  a 
time-reversed  version  of  the  channel  impulse  response,  the  acoustic  propagation  will 
perform  the  matched  filtering  operation,  and  the  received  signal  should  be  temporally 
concentrated.  No  matched  filter  will  be  necessary  at  the  receiver.  They  successfully 
implemented  a  simple  version  of  such  a  system  using  a  reel-to-reel  tape  deck  and 
two  ships  about  18  km  apart  in  2  km  of  water  connected  by  a  radio  link.  The  source 
ship  transmitted  an  impulsive  signal,  and  the  signal  observed  at  the  receiving  ship 
was  recorded  on  the  tape  deck  over  the  radio  channel.  The  source  ship  then  trans¬ 
mitted  a  time-reversed  version  of  the  recorded  impulse  response  by  playing  the  tape 
backwards,  and  the  received  waveform  was  the  desired  temporally-concentrated  sig- 
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nal.  In  addition,  by  making  several  successive  transmissions  using  the  same  reversed 
impulse  response,  Parvulescu  and  Clay  attempted  to  estimate  the  coherence  time 
of  the  channel  from  the  deterioration  of  the  received  impulse.  Unfortunately,  they 
could  not  separate  the  effects  of  channel  variability  from  those  due  to  the  source  ship 
drifting.  The  combination  of  these  effects  resulted  in  the  received  signal  being  almost 
completely  decorrelated  with  the  matched  filter  after  about  30  minutes.  Because  the 
source  ship  had  drifted  approximately  500  m  over  this  time  interval,  Parvulescu  and 
Clay  proposed  this  displacement  was  the  primary  factor  causing  the  mismatch,  but 
acknowledge  without  a  rigorous  control  for  the  experiment  it  would  be  impossible  to 
say  conclusively. 

DiNapoli  et  al.  attempted  a  similar  experiment  over  a  much  longer  range  (50- 
150  km)  and  lower  frequencies  (10-30  Hz)  in  the  Arctic  [48].  Because  the  Arctic 
environment  is  fairly  stable,  they  believed  computational  models  could  predict  the 
impulse  response  of  the  propagation  path,  and  no  attempt  was  made  observe  the 
actual  impulse  response  as  Parvulescu  and  Clay  did.  The  transmissions  DiNapoli  et 
al.  made  using  a  reversed  version  of  the  impulse  response  predicted  by  computational 
models  failed  to  give  a  clearly  identifiable,  temporally  concentrated  arrival.  The  most 
likely  causes  for  this  disappointing  result  are  inaccuracies  in  their  environmental 
knowledge  and  in  the  Arctic  acoustic  models  used  to  predict  the  impulse  response. 
The  impulsive  nature  of  noise  in  the  Arctic  due  to  ice  dynamics  also  would  make  it 
difficult  to  identify  the  desired  arrival  amidst  the  ice  noise. 

The  single  mode  excitation  problem  generalizes  this  approach  from  attempting 
to  create  a  waveform  which  is  concentrated  over  a  short  time  interval  one  which  is 
concentrated  at  a  single  mode.  The  challenge  is  to  determine  which  transmitted 
source  array  weights  will  give  the  desired  mode  coefficients  at  the  feedback  array. 
If  the  channel  is  assumed  to  be  range-invariant,  the  source  array  weights  exciting  a 
single  mode  can  be  computed  from  the  sound  speed  profile  and  bottom  composition 
at  the  source  array.  This  is  analogous  to  assuming  the  channel  is  dispersionless  and 
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has  no  multi-path  interference  for  the  time-domain  experiments  of  Parvulescu  and 
Clay.  Most  shallow  water  channels  are  not  range-invariant.  In  many  range-varying 
channels,  the  weights  computed  assuming  the  channel  is  range-invariant  will  excite  a 
variety  of  modes  at  the  feedback  array  location,  and  not  just  the  desired  single  mode. 
The  source  initially  excites  a  single  mode  which  is  then  distorted  by  propagation 
through  the  channel,  coupling  its  energy  into  several  modes.  The  question  then 
arises  how  best  to  excite  an  initial  distribution  of  modes  at  the  source  array  which  is 
pre-compensated  for  this  distortion  such  that  the  coupling  caused  by  the  channel  will 
cancel  out  all  the  modes  except  the  desired  one  when  the  pressure  field  is  observed 
at  the  feedback  array.  The  feedback  control  algorithm  attempts  to  estimate  the 
matched  signal  for  the  channel  coupling  between  the  source  and  feedback  arrays 
adaptively.  This  multiple  input/multiple  output  problem  is  more  complicated  than 
the  single  source/single  hydrophone  one  addressed  by  Parvulescu  and  Clay,  but  the 
underlying  strategy  of  matching  the  transmitted  signal  to  the  environment  is  the 
same. 


1.6  Summary 

As  noted  at  the  start  of  this  chapter,  the  primary  contribution  of  this  thesis  is  the 
synthesis  of  ideas  from  several  disciplines  to  produce  a  new  algorithm  for  exciting  a 
single  mode.  This  chapter  has  reviewed  the  relevant  background  information  that 
will  be  used  to  derive  the  control  algorithm  in  the  next  chapter.  The  section  on 
mode  filtering  yielded  an  important  perspective  on  the  relationship  among  the  MAP, 
pseudo-inverse,  and  sampled  mode  shape  mode  filters.  The  MAP  filter  was  also 
shown  to  provide  a  theoretical  justification  for  Yang’s  ad  hoc  dropped  eigenvalue 
mode  filter.  Chapter  2  will  show  how  observability  and  reachability  can  provide 
insight  into  how  the  source  and  feedback  array  geometries  affect  the  ability  of  the 
algorithm  to  excite  a  single  mode. 


45 


46 


Chapter  2 


Control  Algorithm 


This  chapter  presents  an  algorithm  for  single  mode  excitation  using  indirect  control 
[14].  An  indirect  control  algorithm  explicitly  estimates  a  model  of  the  plant,  then 
uses  this  model  to  compute  the  control  input  giving  the  desired  behavior.  The 
algorithm  proposed  for  single  mode  excitation  begins  with  some  initial  estimate  of 
the  acoustic  propagation  through  the  feedback  volume.  Based  on  this  estimate,  it 
computes  the  best  set  of  source  array  weights  w  to  get  the  desired  pressure  vector  p^. 
After  exciting  the  channel  using  w,  the  algorithm  compares  the  observed  pressure 
vector  p  to  p,  the  pressure  predicted  by  the  current  channel  estimate.  The  channel 
estimate  is  updated  based  on  the  difference  between  the  observed  and  predicted 
pressure,  and  the  process  begins  again.  Figure  2-1  shows  the  result  of  incorporating 
this  type  of  control  algorithm  into  the  scenario  of  Figure  1-1.  While  this  control 
scheme  is  fairly  simple,  it  clearly  demonstrates  the  potential  of  feedback  control 
for  exciting  a  single  mode.  If  this  straight-forward  scheme  proves  successful,  more 
advanced  control  schemes  like  robust  control  [49]  [50]  [51]  can  be  investigated  for 
this  problem. 

The  control  algorithm  uses  the  replica  matrix  described  in  Section  2.1  to  model 
the  propagation  between  the  source  and  feedback  arrays.  Section  2.1  also  presents  a 
more  precise  definition  of  persistent  excitation  as  first  discussed  in  Section  1.4.  The 
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Figure  2-1:  Control  Algorithm  for  Ocean  Experiment 


source  array  weights  are  computed  from  the  current  estimate  of  the  channel  replica 
matrix  using  the  least  squares  method  described  in  Section  2.2.  The  replica  matrix 
estimator  is  an  important  part  of  the  control  algorithm,  and  this  chapter  presents 
three  different  estimators  for  the  replica  matrix.  Section  2.3  develops  the  condition 
number-limited  least  squares  algorithm  for  estimating  the  channel  replica  matrix. 
The  replica  matrix  estimation  problem  can  also  be  formulated  as  a  Kalman  filtering 
problem,  and  Section  2.4  presents  the  resulting  replica  matrix  estimator.  Section 
2.5  develops  the  least  mean  squares  (LMS)  [13]  estimator  for  the  replica  matrix. 
This  common  adaptive  algorithm  will  be  compared  with  the  other  two  estimation 
algorithms  in  the  experimental  work  of  Chapter  4. 


2.1  Propagation  Model 

The  control  algorithm  presented  in  this  chapter  models  the  pressure  received  at  the 
feedback  array  p  as  a  linear  function  of  the  complex  source  array  weights  w,  i.e., 

P  =  Qw,  (2.1) 

where  Q  is  the  replica  matrix  for  the  channel.  The  column  of  Q  is  the  replica 
vector  at  the  feedback  hydrophone  array  for  the  source  array  element.  Equiva¬ 
lently,  the  (i,  element  of  Q  is  the  Green’s  function  between  the  source  and 
receiver  evaluated  at  the  transmission  frequency.  Considering  the  state  space  mod¬ 
els  of  Eqs.  1.35  and  1.36,  the  replica  matrix  is  Q[n]  =  ’i'rNB[n].  Thus,  this  model 
combines  the  effects  of  the  feedback  array  hydrophone  locations  on  with  the 
transmission  properties  of  B[n]  summarizing  the  source  locations  and  mode  coupling. 
Poor  conditioning  in  Q[n]  can  be  caused  by  any  or  all  of  the  following  phenomena: 
poor  hydrophone  sampling  of  the  modes  at  the  feedback  array,  poor  placement  of 
the  source  array  elements  for  exciting  the  modes  at  the  source  array,  or  unfavorable 
mode  coupling  between  the  source  and  feedback  arrays. 
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The  validity  of  Q[n]  as  a  propagation  model  depends  on  the  time  scales  of  varia¬ 
tions  in  the  channel  relative  to  the  propagation  time  between  the  source  and  feedback 
arrays.  As  noted  in  Section  1.2,  the  replica  matrix  model  is  valid  for  continuous  wave 
(CW)  propagation  in  linear,  time-invariant  propagation  environments.  For  the  typi¬ 
cal  scenario  shown  in  Fig.  1-1,  the  propagation  time  from  source  to  feedback  array  is 
1-2  sec.  The  oceanographic  features  and  processes  dominating  shallow  water  acous¬ 
tics  such  as  internal  waves,  tidal  mixing,  and  bottom  bathymetry  are  believed  to 
vary  on  timescales  much  longer  than  this.  Observations  in  the  recent  SWARM  ex¬ 
periment  indicated  the  predominant  energy  in  the  internal  waves  in  a  typical  coastal 
environment  had  timescales  on  the  order  of  five  to  seven  minutes  [52].  Based  on 
these  observations,  it  appears  safe  to  assume  the  channel  is  time-invariant  over  the 
propagation  interval  between  the  source  and  feedback  arrays. 

The  iteration  time  for  the  algorithm  will  depend  both  on  the  distance  vq  from 
the  source  to  feedback  array  and  the  group  velocities  of  the  modes  as  determined  by 
the  sound  speed  profile  c(r,z).  After  changing  the  source  array  weights,  the  control 
algorithm  must  wait  long  enough  for  the  slowest  mode  to  propagate  from  the  source 
array  to  the  feedback  array  before  observing  the  pressure  field  and  updating  its  esti¬ 
mate  of  the  channel  replica  matrix.  The  control  algorithm  samples  the  propagation 
environment  with  sampling  period  equal  to  its  iteration  time.  The  Nyquist  sampling 
theorem  [53]  indicates  that  setting  the  iteration  time  to  be  at  least  the  propagation 
time  for  the  slowest  mode  means  that  oceanographic  processes  at  timescales  shorter 
than  twice  this  propagation  time  will  be  aliased  in  the  observations  of  the  environ¬ 
ment.  This  best  case  limit  on  the  bandwidth  of  observable  processes  must  be  taken 
into  account  when  designing  the  experiment.  Due  to  the  averaging  and  windowing 
of  data  employed  by  the  estimators  described  in  this  chapter,  the  limit  on  the  fastest 
processes  the  control  algorithm  can  expect  to  track  is  much  slower  than  the  limit 
specified  by  the  Nyquist  sampling  theorem. 

A  sequence  of  inputs  w[0], . . . ,  w[n]  is  said  to  be  persistent  if  it  allows  unique 
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identification  of  one  set  of  model  parameters  out  of  all  possible  models  in  a  given 
model  class  [44].  Since  the  replica  matrix  model  for  propagation  defines  the  space  of 
possible  models  to  be  all  AT  x  L  complex  matrices  persistency  can  be  defined 

more  precisely  than  the  general  discussion  in  Section  1.4.  To  identify  a  unique  Q, 
the  sequence  of  input  vectors  must  span  C^.  Thus,  persistency  over  the  interval  [0,  n] 
is  equivalent  to  [w[0]|w[l]| . . .  |w[n]]  having  rank  L.  This  condition  is  also  equivalent 

n 

to  the  square  matrix  ^w[z]w’®[i]  having  full  rank. 

t=0 

Because  the  control  algorithm  attempts  to  hold  the  pressure  field  observed  at 
the  feedback  array  fixed,  the  input  signal  usually  will  not  vary  after  an  initial  con¬ 
vergence  transient  if  the  replica  matrix  estimate  Q  is  accurate  and  the  channel  is 
quasi-stationary.  Later  discussions  will  illustrate  how  this  impersistent  or  nearly 
impersistent  input  sequence  affects  the  convergence  and  transient  behavior  of  the 
estimators.  If  the  steady  state  source  array  weights  are  orthogonal  to  the  er¬ 
rors  in  the  rows  of  Q[n],  these  errors  will  not  be  revealed  by  the  algorithm  in  the 
steady  state.  Although  the  errors  in  Q  are  undetectable  by  w„s,  may  not  be  the 
source  array  weight  vector  giving  the  best  performance.  The  control  algorithm  does 
have  fixed  points  other  than  the  global  optimum  w/,5  =  Q^p^,  and  thus  may  con¬ 
verge  to  a  non-optimal  solution  for  w.  In  the  simulations  and  experiments  described 
later  in  this  thesis,  the  performance  of  the  algorithm  was  limited  by  the  observation 
noise,  and  not  by  convergence  to  non-optimal  solutions.  The  fixed  points  at  these 
non-optimal  solutions  are  generally  unstable  equilibria,  so  observation  noise  should 
provide  sufficient  impetus  to  propel  the  control  algorithm  away  from  them.  This 
likely  explains  why  no  evidence  was  observed  of  the  algorithm  getting  stuck  at  non¬ 
optimum  source  weights  in  the  simulations  or  laboratory  waveguide  experiments. 


51 


2.2  Source  Array  Weight  Selection 

The  least  squares  solution  to  the  system  of  linear  equations 


Pd  =  Qw 


(2.2) 


yields  the  best  source  array  weight  vector  for  obtaining  the  desired  pressure  Pd  given 
the  current  replica  matrix  Q.  Assuming  N  >  L,  the  singular  value  decomposition  of 
Q  is 


Q  =  Uq 


" 

O'Ql 

0 

0 

0 

TT'  H 

VQ  > 


(2.3) 


where  Uq  is  iV  x  AT,  Vq  is  L  x  L,  the  partitioned  matrix  on  the  right-hand 
side  is  N  X  L,  and  the  diagonal  elements  of  that  matrix  are  ordered  such  that 
>  o'Q2  >  •  ■  •  >  o'QL  >  0  [37].  The  pseudo-inverse  of  Q,  is  defined  to  be 


Q^  =  Vq, 


- 

^Q1 

0 

0 

0 

u, 


H 


(2.4) 


except  that  any  crgi’s  equal  to  zero  are  not  inverted,  but  left  as  zero.  If  all  the 
singular  values  aqj  are  non-zero,  w  =  Q^Pd  gives  the  unique  source  array  weight 
vector  minimizing  the  error  between  the  desired  and  obtained  pressures.  If  the 
propagation  is  completely  described  by  M  trapped  modes,  and  M  <  L,  there  will  be 
singular  values  of  Q  that  are  zero,  and  the  solution  minimizing  the  magnitude  of  the 
error  is  not  unique.  In  this  case,  w  =  Q^Pd  gives  the  minimum  norm  source  array 
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weight  vector  with  the  minimum  squared  error  between  the  observed  pressure  p  and 
the  desired  pressure  for  the  true  channel  replica  matrix  Q.  The  weight  selection 
step  uses  the  same  method  to  find  w  except  Q  is  replaced  by  the  current  replica 
matrix  estimate  Q.  For  the  iteration,  the  source  array  weights  are  determined 
by 


w[n]  =  (Q[n])^ 


Prf) 


(2.5) 


where  Q[n]  is  the  estimate  of  the  channel  replica  matrix  at  the  iteration.  This 
gives  the  source  array  weight  vector  that  would  be  the  best  fit  if  the  current  channel 
estimate  QM  were  correct. 

The  weight  selection  step  proposed  in  Eq.  2.5  minimizes  the  error  in  the  pressure 
domain,  which  does  not  generally  give  the  minimum  error  in  the  mode  domain.  It  is 
possible  to  compute  w  to  minimize  the  error  in  the  mode  domain  using  an  estimate 
of  B[n]  from  the  state-space  model  of  Eq.  1.35.  The  weight  selection  step  then  uses 
solves  Eq.  1.35  in  the  least-squares  sense  assuming  the  state  noise  f[n]  is  zero.  The 
resulting  source  array  weight  vector 


w[n]  =  (B(n|)*dj, 


(2.6) 


should  minimize  the  error  in  the  mode  domain.  In  practice,  this  approach  can  be 
problematic  for  numerical  reasons  discussed  below. 

It  is  necessary  to  develop  a  precise  framework  to  discuss  the  important  issues 
that  arise  when  some  of  the  matrices  in  the  control  algorithm  computations  are  not 
mathematically  singular,  but  have  very  large  condition  numbers.  Golub  and  Van 
Loan  define  the  e-rank,  r^,  of  the  matrix  Q  to  be  the  number  of  singular  values  oi  >  e, 
i.e.,  if  re  ==  rank(Q,  e)  =  m,  then  oqi  >  (7q2  >  . . .  >  (Jq^  >  e  >  aQ,„+i  >  gqi  >  0 
[37].  It  is  straight-forward  to  extend  this  idea  and  define  the  e-span  of  a  matrix  with 
e-rank  Tj  as  the  span  of  {ui, . . . ,  Ur^},  where  the  Uj  are  the  columns  of  Uq,  the  left 
orthogonal  matrix  of  the  SVD.  Similarly,  the  e-nullspace  of  the  matrix  can  be  defined 
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to  be  the  space  spanned  by  {vr^+i, . . . ,  v/,},  where  the  Vj’s  are  the  columns  of  Vq, 
the  right  orthogonal  matrix  of  the  SVD. 

If  Q[n]  has  e-rank  less  than  M  for  some  e  appropriate  to  the  scale  of  the  exper¬ 
iments,  it  indicates  the  system  is  either  numerically  unreachable  or  unobservable, 
either  of  which  may  be  problematic  for  mode  excitation.  The  terms  nearly  unreach¬ 
able  and  nearly  unobservable  will  be  used  to  describe  these  numerically  troublesome 
conditions.  If  the  system  is  nearly  unreachable,  both  the  pressure  and  mode  domain 
approaches  for  computing  w[n]  may  have  difficulty  exciting  the  desired  mode.  In 
practical  terms,  this  condition  corresponds  to  the  scenario  where  the  source  element 
locations  and  propagation  conditions  in  the  feedback  volume  combine  unfavorably 
so  that  some  modes  cannot  be  excited  to  any  significant  level  given  the  power  lim¬ 
itations  on  the  sources.  If  the  desired  pressure  e  e-span{Q},  then  the  desired 
pressure  field  may  be  obtained  within  the  power  limits  of  the  practical  system  as 
refiected  by  the  choice  of  e.  Otherwise,  the  power  required  to  excite  this  particular 
mode  exceeds  the  limitations  of  the  sources.  It  is  possible  that  for  a  given  source 
array  geometry  and  propagation  conditions  which  result  in  numerical  unreachability, 
some  modes  will  still  be  excitable  and  others  will  not,  depending  on  the  e-span  of  Q. 

If  the  e-rank  deficiency  of  Q[n]  is  because  the  system  is  nearly  unobservable, 
there  may  be  undesired  trapped  modes  propagating  at  the  feedback  array  even  when 
the  observed  pressure  p[n]  is  a  close  fit  to  the  desired  pressure  p^.  Numerical  un¬ 
observability  can  be  due  to  one  or  more  causes.  Some  modes  may  be  sampled  at 
or  close  to  a  pressure  null,  so  that  any  energy  in  these  modes  cannot  be  seen  over 
the  observation  noise.  Alternatively,  if  one  the  columns  of  the  sampled  mode  shapes 
can  be  written  as  a  linear  combination  of  other  columns,  the  pressure  field  cannot 
be  uniquely  resolved  into  modes.  These  conditions  are  equivalent  to  the  poorly  con¬ 
ditioned  sampling  for  ^  discussed  in  Section  1.3.  If  the  feedback  array  sufficiently 
samples  the  water  column  snch  that  the  system  is  observable,  the  issue  is  not  as  im¬ 
portant.  In  this  case,  exciting  a  pressure  field  close  to  the  desired  p^  will  correspond 
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to  a  mode  coefficient  vector  d[n]  close  to  the  desired  d^. 

The  weight  vector  minimizing  Ijpa— p[n]||  will  not  generally  minimize  ||d(i— d[n]||. 
The  errors  in  the  two  domains  are  related  by  a  linear  but  generally  non-unitary 
transformation 

Pd-PN  =  ^(dd-d[n]).  (2.7) 

To  see  why  the  minimum  error  in  the  pressure  domain  does  not  necessarily  correspond 
to  the  minimum  error  in  the  mode  domain,  consider  the  extreme  case  when  ^ 
possesses  a  non- trivial  null-space,  i.e.,  when  the  system  is  unobservable.  Let  d[n] 
be  the  mode  vector  minimizing  the  error  in  the  mode  domain,  and  p[n]  be  resulting 
pressure  vector,  i.e.,  p[n]  =  ’®'[n]d[n].  Any  d[n]  =  d[n]  -f-  Ad,  where  Ad  is  in  the 
null-space  of  ’i'  will  attain  the  same  pressure  domain  error,  but  not  minimize  the 
mode  domain  error.  The  argument  is  somewhat  more  involved  when  ’i'  does  not 
have  a  null-space,  but  so  long  as  the  singular  values  of  'i'  are  not  all  equal,  the 
problem  exists.  If  the  system  is  well-sampled  such  that  the  o-j’s  of  are  at  least 
commensurate  with  each  other  if  not  equal,  the  weight  vector  giving  the  minimum 
error  in  the  pressure  domain  will  give  a  small  if  not  minimum  error  in  the  mode 
domain.  In  the  limiting  case  with  the  array  spanning  the  entire  water  column  as 
N  oo,  the  CTsf  t’s  approach  one  and  the  solution  minimizing  the  error  in  one  domain 
also  minimizes  it  in  the  other.  As  noted  in  Section  1.3  the  difficulty  of  observing  the 
pressure  in  the  bottom  prevents  the  singular  values  from  actually  reaching  one,  but 
generally  they  will  come  very  close  to  it. 

To  see  more  clearly  the  difficulty  in  explicitly  estimating  d[n]  when  ’i'  is  nearly 
singular,  consider  the  scenario  where  the  observation  noise  n  has  a  Gaussian,  spa¬ 
tially  white  component  with  variance  cr^,  and  the  propagation  physics  ensure  the 
magnitudes  of  all  mode  coefficients  are  bounded  by  some  constant  dmax,  i-e., 

niax||dj[n]||  <  (2.8) 

t,n 
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If  rank(’®',  £r„/Arf)  ^  min(M,iV),  the  observations  d[n]  =  will  be  dominated 

by  noise.  One  or  more  of  the  on  the  diagonal  of  will  be  large  enough  that  it 
will  make  the  white  noise  component  in  the  Vj  direction  larger  than  the  coefficients 
of  the  propagating  modes.  This  will  slow  the  convergence  of  the  algorithm  since 
the  estimator  will  have  to  average  many  observations  of  d[n]  to  obtain  an  accurate 
estimate  of  B  from  these  low  signal-to-noise  ratio  observations.  Because  B[n]  will 
converge  slowly  and  fluctuate  significantly  due  to  noise,  the  weight  vector  given  by 
Eq.  2.6  will  vary  considerably.  Under  these  conditions,  the  pressure  vector  observed 
at  the  hydrophone  will  be  unsatisfactory  for  two  reasons.  Because  the  system  is 
nearly  unobservable,  undesired  modes  may  be  propagating  at  the  feedback  array 
which  are  not  visible  over  the  observation  noise.  In  addition,  fluctuations  in  B[n]  will 
make  w[n]  and  consequently  p[n]  vary  significantly  from  one  iteration  to  the  next. 
The  pressure  field  at  the  feedback  array  may  not  only  contain  undesired  modes,  but 
may  not  even  consist  of  the  same  modes  on  successive  iterations.  For  these  reasons, 
the  pressure  at  the  feedback  array  will  be  a  poor  source  signal  for  oceanographic 
observers  further  downrange  at  the  observation  array. 

If  the  modes  are  nearly  unobservable,  formulating  the  problem  in  the  pressure 
domain  will  not  alter  the  fact  that  the  pressure  field  may  consist  of  many  modes, 
rather  than  the  desired  single  mode.  If  the  observation  noise  n[n]  is  relatively  weak 
compared  to  the  pressure  due  to  the  propagating  modes,  the  pressure  observed  at 
the  feedback  array  will  not  vary  much  when  w[n]  stays  fixed.  Because  the  noise  is 
not  amplified  as  it  was  when  estimating  d[n],  the  observations  p[n]  are  more  reliable, 
and  the  estimate  Q[n]  will  converge  faster  than  B[n]  and  be  less  sensitive  to  noise. 
In  this  case  the  pressure  domain  approach  will  give  a  pressure  field  at  the  feedback 
array  which  is  at  least  consistent  over  time  even  if  it  does  not  give  the  desired  single 
mode.  Although  from  the  perspective  of  the  oceanographic  observers  a  high-fidelity 
single  mode  is  the  best  source,  a  consistent  mode  distribution  of  lower  fidelity  is 
preferable  to  the  erratic  mode  coefficients  obtained  by  computing  w[n]  from  B[n] 
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when  ’i'[n]  is  poorly-conditioned. 

If  the  spread  of  singular  values  of  ^  is  such  that  the  modes  are  practically  ob¬ 
servable,  i.e.,  (rank('y?,  cr„/Ad)  =  M),  but  the  -wfn]  minimizing  ||pd  -  p[n]j|  gives 
unacceptably  large  errors  in  the  mode  domain,  the  weight  selection  step  can  be 
reformulated  as  a  weighted  least  squares  problem  using  an  error  metric  where  all 
components  in  the  mode  domain  are  equally  weighted.  The  source  array  weights  w 
are  chosen  to  minimize 

||e[n]|p  =  ||Re[n]|p  =  (p[n]  -  pd)“R2(p[n]  -  pd),  (2.9) 

where  R  is  an  AT  x  AT  positive  semi-definite  Hermitian  matrix.  Choosing 


(2.10) 


gives  the  weighted  pressure  domain  error  whose  minimization  corresponds  to  the 
minimum  error  in  the  mode  domain.  Note  that  R  has  the  same  orthogonal  comple¬ 
ment  as  indicating  pressure  domain  errors  orthogonal  to  the  columns  of  are 
irrelevant  to  this  metric.  This  is  reasonable  since  if  ^  is  accurate,  errors  in  these 
directions  could  not  be  due  to  the  modes  and  must  be  observational  noise.  Solving 
the  resulting  least  squares  problem  gives  the  new  criteria  for  selecting  the  source 
array  weight  vector: 


R=U^ 


.-1 

'4-1 


0 


0 


-1 


0 


0 


0 


w[n]  =  (Q“[n]R^Q[n])^  Q[n]R^Pd. 


(2.11) 
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If  the  replica  matrix  estimate  Q[7i]  =  Eq.  2.11  reduces  to 

w(nj  =  (B[n])'  (s'pj)  =  (B[ti])*  dj  (2.12) 

which  would  be  the  weight  vector  minimizing  the  error  in  the  mode  domain  if  an  ex¬ 
plicit  estimate  of  B[n]  were  available  to  the  weight  selection  step.  The  weight  vector 
that  minimizes  the  weighted  pressure  domain  error  defined  with  Eqs.  2.9  and  2.10 
gives  the  desired  weight  vector  minimizing  the  error  in  the  mode  domain. 

This  section  has  developed  a  method  for  estimating  the  best  source  array  weight 
vector  w[n]  to  excite  the  desired  pressure  profile.  This  estimate  depends  on  the 
current  estimate  of  the  channel  replica  matrix  Q[n].  The  remaining  three  sections 
of  the  chapter  develop  three  different  methods  for  estimating  Q[n]  for  the  feedback 
volume. 


2.3  Condition  Number-limited  Least  Squares 

The  source  array  weight  vector  w[n]  is  a  function  of  the  current  estimate  of  the  replica 
matrix,  Q[n].  In  order  to  determine  w[n],  the  control  algorithm  must  compute  Q[n], 
the  estimate  of  Q  at  the  iteration,  from  the  information  available  about  the 
channel,  specifically,  the  sequence  of  source  array  weights  w[0], . . . ,  w[n  —  1]  and  the 
pressures  observed  at  the  feedback  array  p[0], . . . ,  p[n  —  1].  Each  of  the  N  rows  of 


qiN 


Q[n]  = 


(2.13) 


qvW 


can  be  estimated  independently  of  the  other  rows  using  the  observed  pressure  at  the 
hydrophone  of  the  feedback  array  corresponding  to  that  row  of  the  replica  matrix. 
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If  the  rows  of  Q[n]  are  assumed  to  be  constant,  the  sequence  of  source  array 
weights  and  feedback  array  pressures  form  a  system  of  linear  equations: 

bt[0]| . . .  \pi[n  -  1]]  =  Qi  [w[0]| . . .  |w[n  -  1]] ,  (2.14) 

or 

Pi[n-l]  =  q,W[n-l],  (2.15) 

where 

Pi[n-1]  =  bi[0]l---bib- 1]]  (2-16) 

Wb  —  1]  =  [w[0]| . . .  |w[n  —  1]] .  (2-17) 

The  least  squares  solution  to  this  system  provides  q^n] 

qi[n]  =  Pi[n  -  l]W“b  -  1]  (w[n  -  l]W“[n  -  1])~^ .  (2.18) 

This  equation  can  be  manipulated  into  a  recursive  form  by  following  the  devel¬ 
opment  of  the  Recursive  Least  Squares  (RLS)  algorithm  [54].  Although  the  matrix 
inversion  lemma  used  in  that  algorithm  is  a  tempting  means  of  avoiding  the  inversion 
in  Eq.  2.18,  the  lemma  is  too  numerically  sensitive  to  be  useful  for  this  problem  in 
practice.  The  extreme  conditioning  of  W[n  —  l]W^[n  —  1]  and  the  numerical  issues 
associated  with  it  will  be  discussed  in  more  detail  later  in  this  section.  Because 
W[n  —  l]W**[n  —  1]  is  not  a  function  of  i,  the  row  number,  the  matrix  inverted  is 
the  same  for  each  row  of  Q.  This  allows  all  rows  of  the  replica  matrix  estimate  to 
be  computed  simultaneously.  The  same  algebraic  manipulations  used  to  derive  the 
RLS  algorithm  in  [54]  yield  a  partially  recursive  form  of  the  least  squares  estimator 
for  the  channel  replica  matrix: 
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n  +  1] 

=  Q[n]  -1-  a[n]k”[n] 

(2.19) 

a[n] 

=  pN  -  QNw[n] 

(2.20) 

k[n] 

=  #“^[n]w[n] 

(2.21) 

s 

e 

+ 

1 

II 

II 

(2.22) 

Although  this  estimator  must  invert  $[n]  for  each  iteration,  it  is  recursive  in  the 
sense  that  computing  Q[n  +  1]  requires  only  the  current  observations  (p[n]  and 
w[n]),  the  previous  estimate  Q[n]  and  the  source  weight  correlation  matrix  $[n], 
itself  easily  updated.  The  complete  sequence  of  source  array  weights  and  feedback 
array  pressures  need  not  be  saved.  Note  that  a[n]  is  the  error  signal  between  the 
observed  pressure  and  the  pressure  predicted  by  the  current  estimate  Q[n]. 

The  algorithm  must  be  initialized  with  a  replica  matrix  Q[0]  and  source  weight 
correlation  #[0].  These  can  be  obtained  by  sequentially  turning  on  each  element  of 
the  source  array  with  gain  1,  then  collecting  the  observed  pressures  into  the  columns 
of  Q[0]  and  setting  $[0]  =  I.  Alternatively,  Q[0]  can  be  set  equal  to  some  a  priori 
guess  at  the  replica  matrix,  and  ^[0]  =  for  some  small  constant  5. 

In  realistic  ocean  environments,  Q  will  not  usually  be  constant,  but  slowly  vary¬ 
ing.  The  estimator  specified  in  Eqs.  2.19-2.22  can  be  modified  to  include  an  expo¬ 
nential  forgetting  factor  7  to  depreciate  older  observations  and  allow  the  estimator 
to  respond  to  changes  in  the  ocean  environment  [54].  This  forgetting  factor  can  be 
incorporated  into  the  estimator  by  modifying  Eq.  2.22  to 

$[n]  =  ')^[n  —  1]  +  w[n]w**[n],  (2.23) 

where  0  <  7  <  1.  Multiplying  by  7  at  each  iteration  modifies  the  estimator  to 

n 

choose  qt[n]  to  minimize  5^7”''^||Pi[^]  —  This  new  error  criteria  reduces 

1=0 
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the  importance  of  old  observations  until  they  are  eventually  insignificant.  The  best 
choice  of  7  will  depend  on  the  bandwidth  and  intensity  of  the  oceanographic  pro¬ 
cesses  affecting  acoustic  propagation  and  the  iteration  time  of  the  control  algorithm. 
Choosing  7  to  be  too  small  will  not  allow  the  estimator  to  average  enough  observa¬ 
tions  to  reduce  the  effects  of  observation  noise.  Too  large  a  choice  for  7  will  not  allow 
the  estimator  to  respond  expeditiously  to  changes  in  the  propagation  conditions.  The 
choice  of  7  must  compromise  between  these  conflicting  demands. 

The  discussion  of  persistency  in  Section  2.1  observed  that  single  mode  excitation 
will  generally  require  w[n]  to  be  constrained  within  a  small  region  of  C^.  The  con¬ 
dition  for  a  persistent  input  sequence  given  there  is  equivalent  to  ^[n]  as  defined 
in  Eq.  2.22  having  full  rank.  An  impersistent  input  causes  problems  for  this  algo¬ 
rithm  because  as  w[n]  remains  constant  or  nearly  constant,  #[n]  will  become  poorly 
conditioned.  This  causes  difficulty  computing  before  #[n]  is  mathematically 

singular. 

The  large  condition  number  of  #[n]  will  cause  the  estimator  to  become  overly 
sensitive  to  observation  noise  even  before  problems  computing  the  inverse  arise.  A 
large  sample  of  noise  may  make  the  error  vector  a[n]  7^  0  even  when  Q[n]  is  accurate. 
The  resulting  small  perturbation  in  Q[n]  will  cause  a  slight  perturbation  in  w[n] 
away  from  the  correct  source  array  weight  vector.  As  the  subsidiary  eigenvalues  of 
$[n]  grow  arbitrarily  small,  even  a  minor  perturbation  in  w[n]  from  its  fixed  value 
will  cause  k[n]  to  have  a  very  large  norm.  This  could  result  in  a  large  correction 
to  Q[n],  and  subsequently  a  large  error  in  w[n]  and  p[n]  for  the  next  iteration. 
Although  the  estimator  will  quickly  recover  from  this  noise-induced  transient,  these 
excursions  produce  undesirable  spikes  in  the  error  between  the  desired  and  observed 
mode  coefficients.  These  spikes  detract  from  the  feedback  array’s  value  as  a  single 
mode  source  to  observers  further  downrange  in  the  scenario  of  Figure  1-1. 

These  transients  may  be  prevented  by  limiting  the  condition  nnmber  of  the  inverse 
of  ^[n]  to  be  no  greater  than  some  threshold  77.  To  define  this  condition  number- 
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limited  inversion  more  precisely,  let  the  eigenvalue  decomposition  of  #[n]  be 


^[n]  = 


0 


0 


(2.24) 


where  A$i  >  A^2  >  •  •  •  >  >  0.  The  condition  number-limited  inverse  ^cz,W 

defined  to  be 


V^«,  (2.25) 


where  Aj  =  max(Ai,  Ai/77).  This  insures  the  condition  number  of  never 

greater  than  rj,  while  Ai,  the  eigenvalue  corresponding  to  the  predominant  eigen¬ 
vector,  is  correctly  inverted.  Using  ^c\X^  compute  k[n]  insures  that  no  compo¬ 
nent  of  w[n]  is  unduly  amplified,  and  thus  the  norm  of  k[n]  never  grows  too  large. 
Replacing  by  Eq.  2.21  completes  the  definition  of  the  condition 

number-limited  least  squares  (CNLLS)  algorithm.  Note  that  the  source  weight  cor¬ 
relation  matrix  4^[n]  propagates  with  unlimited  condition  number.  If  the  channel  or 
desired  mode  changes  abruptly,  the  transient  response  is  faster  if  $[n]  has  not  been 
condition  number-limited.  Appendix  A  contains  a  simple  example  demonstrating 
the  importance  of  limiting  the  condition- number  of  the  inverse  and  not  $[n]  itself. 

Another  consequence  of  an  impersistent  excitation  is  that  the  CNLLS  estimate  of 
the  replica  matrix  may  converge  to  an  inaccurate  estimate.  From  Eqs.  2.19-2.22,  it 
is  apparent  that  any  value  of  Q[n]  yielding  a[n]  =  0  is  a  fixed  point  of  the  estimator. 
Writing  p[n]  =  Q[n]w[n]  reveals  that  a[n]  =  0  is  equivalent  to  (Q[n]  —  Q[n])w[n]  =  0 
or  (AQ[n])w[n]  =  0  where  AQ[n]  is  the  error  in  the  current  replica  matrix  estimate. 
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This  condition  does  not  imply  AQ[n]  =  0,  but  only  that  the  error  in  each  row 
of  the  estimate  is  orthogonal  to  w[n].  Consequently,  (Q[n]  -  AQ[n])lprf  does  not 
necessarily  give  the  w[n]  minimizing  ||p[n]  -  pd||.  As  discussed  in  Sec.  2.1,  these 
fixed  points  are  unstable  equilibria,  and  in  practice  the  observation  noise  appears 
to  prevent  the  algorithm  from  stalling  at  these  incorrect  channel  replica  estimates. 
The  resulting  jitter  in  w[n]  keeps  the  input  persistent  enough  to  prevent  any  severe 
errors  in  QM- 

2.4  Kalman  Filter  Estimator 

The  rows  of  Q[n]  can  also  be  estimated  as  the  state  vectors  for  N  separate  Kalman 
filters.  In  order  to  apply  the  Kalman  filter  to  the  problem,  the  rows  of  Q[n]  are  mod¬ 
eled  as  evolving  as  separate  first-order  Gauss-Markov  processes  [55].  The  covariance 
of  the  Gaussian  process  driving  the  state  equations  must  be  observed,  computed  from 
oceanographic  models,  or  assumed.  If  the  rows  of  Q[n]  propagate  independently,  the 
state  space  equations  can  be  written  as 

qi[n  +  1]  =  aiCii[n\  +  fj[n]  (2.26) 

Pi[n]  =  qt[n]w[n] -fnjn],  (2.27) 

where  fjfn]  is  the  Gaussian  random  process  driving  the  state  update  equations  with 
covariance  Pfjn],  pi[n\  is  the  pressure  observed  at  the  hydrophone  of  the  feedback 
array,  and  ni[n]  is  the  observation  noise  at  that  hydrophone,  with  variance  . 
Note  the  assumption  of  independence  is  very  conservative  in  this  context,  as  the 
estimators  do  not  attempt  to  exploit  any  correlation  between  the  rows  of  Q.  This 
also  prevents  the  estimators  from  suffering  any  deterioration  in  performance  due  to 
possible  mismatch  between  the  true  values  of  the  correlations  and  the  values  used 
by  the  estimator. 

Equations  2.26  and  2.27  are  identical  to  the  formulations  of  the  state-space  model 
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in  [43]  and  [55],  except  the  format  has  been  appropriately  modified  to  make  qt[n] 
a  row  vector  instead  of  a  column  vector.  Exploiting  the  development  of  [43],  the 
Kalman  filter  estimate  of  qt[»^]  can  be  written  as 
Prediction  Step 


qi[n]n  -  1]  =  a;iqt[n  -  l]n  -  1]  (2.28) 

P[n]n  —  1]  =  a!iP[n  —  l]n  —  1]  +  Pf;  (2.29) 


Update  Step 

(2.30) 

(2.31) 

where  qi[n]n]  and  qf[n]n— 1]  are  E  {q[n])pi[0], . .  ■  ,Pi[n]}  and  E  {q[n]|pi[0], . .  -  ,Pt[’^  —  1]} 
respectively.  The  covariance  matrices  are  defined  as 


qi[n|n] 

P[n]n] 


qi[nln  -  1]  + 

P[n]n  —  1]  — 


(pi[n]  -  qi[n]n  -  l]w[n])w^[n]P[nln  -  1] 
w^[n]P[n]n  —  l]w[n]  +  cr^.  ’ 

P[n]n  —  l]w[n]w'^[n]P[nl?2  —  1] 
w^[n]P[n]n  —  l]w[n]  + 


P[n]n]  =  cov  {q[n]  —  q[n]n]} 
P[n]n— 1]  =  cov  {q[n]  —  q[n]n  —  1]}  . 


(2.32) 

(2.33) 


If  Pf^  and  Oj  are  assumed  to  be  identical  for  all  the  rows  of  Q,  and  the  observation 
noise  n[n]  is  spatially  white  and  Gaussian  with  variance  cr^,  the  row  estimates  may 
all  be  computed  simultaneously  as 
Prediction  Step 


Q[n]n  —  1]  =  Q:Q[n  — Ijn— 1]  (2-34) 

P[nln-1]  =  o2p[n- l]n- 1] +Pf  (2.35) 
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Update  Step 


Q[n|n] 

P[nln] 


Q[n|n  -  1]  + 
P[n|n  —  1]  — 


(pN  -  Q[n|n  -  l]w[n])w‘^[n]P[n|n  -  1] 
wH[n]P[n|n  -  l]w[n]  +  cr^ 

P[n|n  —  l]w[n]w^[n]P[n|n  —  1] 
w*^[n]P[n|n  —  l]w[n]  + 


(2.36) 

(2.37) 


The  Kalman  filter  algorithm  begins  with  initial  estimates  Q[0|0]  and  P[0|0].  These 
may  be  obtained  as  specified  for  the  CNLLS  algorithm  by  sequentially  activating  each 
element  of  the  source  array,  or  estimated  from  a  priori  knowledge  of  the  environment. 

The  Kalman  filter  estimator  has  an  advantage  over  the  CNLLS  algorithm  because 
the  former  provides  mechanisms  for  incorporating  and  exploiting  any  additional  in¬ 
formation  available  about  the  spatial  correlation  of  the  noise  or  replica  matrix  evolu¬ 
tion  over  time.  For  instance,  if  the  observation  noise  is  dominated  by  sea  surface  noise 
matching  the  Kuperman-Ingenito  model,  the  cr^.  terms  can  be  set  appropriately  for 
each  row.  Currently,  no  detailed  and  experimentally- verified  models  exist  describing 
how  the  statistics  of  replica  matrices  are  affected  by  oceanographic  processes  such  as 
internal  waves  and  tidal  mixing.  If  experiments  using  single  mode  propagation  and 
other  techniques  yield  such  models,  this  information  can  be  incorporated  into  the 
estimator  in  two  ways.  First,  these  models  could  give  predictions  for  Pf..  Second,  in 
many  conditions  such  as  internal  wave  propagation  through  the  control  feedback  vol¬ 
ume,  the  updating  of  the  Qj’s  will  not  be  independent.  An  accurate  model  of  internal 
wave  dynamics  could  estimate  the  correlation  between  the  rows  of  Q[n].  If  the  rows 
are  combined  into  one  large  1  x  NL  state  vector,  these  correlation  estimates  would 
give  Pf  for  this  new  state  vector,  and  the  Kalman  filter  would  exploit  the  correla¬ 
tion  structure  of  these  rows  in  updating  its  estimate  of  Q[n].  The  Kalman  filter’s 
ability  to  assimilate  such  information  as  it  becomes  available  recommends  it  over 
the  least  squares  approach.  As  noted  earlier,  the  possible  benefits  of  including  such 
information  in  the  model  comes  with  the  risk  of  concomitant  mismatch  problems  if 
the  information  is  incorrect. 
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Two  recent  papers  have  presented  preliminary  results  on  the  statistics  and  dy¬ 
namics  of  interval  wave  processes.  Lynch  et  al.  [56]  mecisured  travel-time  variation 
in  acoustic  propagation  due  to  large  scale  internal  tides  in  the  Barrents  Sea  Frontal 
Region.  Unfortunately,  equipment  failure  prevented  the  determination  of  directional 
spectra  of  the  internal  wave  field,  and  also  scuttled  plans  to  observe  propagation 
along  a  path  which  would  have  isolated  internal  wave  effects  on  acoustics  from  the 
effects  of  frontal  dynamics  and  bottom  layer  fluid  dynamics.  Candy  and  Chambers 
proposed  a  discrete-time  state-space  model  for  internal  wave  dynamics  [57].  They 
used  this  model  to  enhance  the  signature  of  internal  waves  observed  by  a  horizon¬ 
tal  current  meter  array.  The  state  space  model  required  careful  hand-tuning  of  the 
covariance  matrices  for  the  specific  environment  and  waves  in  order  to  detect  the 
internal  waves.  As  more  data  is  obtained  about  internal  wave  dynamics,  a  careful 
evaluation  of  the  potential  benefits  versus  risks  must  be  made  before  incorporating 
the  model  statistics  into  the  Kalman  filter. 

The  Kalman  filter  estimator  suffers  when  the  input  sequence  is  not  persistent. 
The  input  sequence  w[n]  plays  the  role  of  the  output  matrix  C[n]  in  Eq.  1.34.  A 
system  is  defined  to  be  uniformly  completely  observable  if  there  exist  a.  and  m  such 
that  the  smallest  eigenvalue  of  the  observability  gramian  M[n,  n  +  m]  is  greater 
than  a  for  all  n  [43].  The  Kalman  filter  to  estimate  the  state  of  a  system  in  the 
form  of  Eqs.  1.33  and  1.34  will  have  exponentially  stable  dynamics  if  the  pair  of 
matrices  ^A[n],  A'“,[/^[n]C[n]j  are  uniformly  completely  observable,  and  the  pair 
(A[n],  Kfr[n])  are  uniformly  completely  observable.  The  former  criteria  is  the  crit¬ 
ical  one  for  the  Kalman  filter  replica  matrix  estimator.  For  the  state  space  model 
described  in  Eqs.  2.26  and  2.27,  the  observability  criteria  applies  to  (ajl,  cr^/w[n]  j , 
giving  an  observability  gramian 

n-{-m 

M[n, n  +  m]  =  (2.38) 

i=n 


66 


If  w[n]  converges  to  a  fixed  value,  this  matrix  will  not  satisfy  the  condition  for 
uniform  observability.  If  |at|  <  1,  it  will  only  exacerbate  this  conditioning  problem 
since  w[n]’s  earlier  excursions  away  from  the  final  value  will  be  weighted  even  less. 
The  control  algorithm’s  tendency  to  produce  an  impersistent  input  sequence  means 
the  Kalman  filter  cannot  theoretically  be  guaranteed  to  converge.  However,  the 
algorithm  had  no  difficulty  converging  in  either  simulations  or  laboratory  waveguide 
tests. 


2.5  Least  Mean  Squares  Estimator 

The  Least  Mean  Square  (LMS)  Algorithm  [13],  [54]  is  a  computationally  simple  and 
commonly  used  algorithm.  In  terms  of  the  variables  of  the  replica  matrix  estimation 
problem,  the  algorithm  adaptively  chooses  a  row  vector  q[n]  to  give  the  best  least 
squares  solution  to 

Pi[n]  =  qi[n]w[n].  (2.39) 

The  LMS  algorithm  approximates  the  steepest  descent  algorithm  by  using  the  instan¬ 
taneous  values  ofpt[n]  and  w[n]  to  estimate  the  cross-correlation  vector  E  |pi[n]w®[n] 
and  autocorrelation  matrix  E  |w[n]w^[n]|.  The  resulting  estimator  is 

qi[n  -I- 1]  =  qi[n]  4-  n  {pi[n]  -  qi[n]'w[n])  w’^[n],  (2.40) 

where  /x  is  a  scalar  controlling  the  correction  step  size  at  each  iteration.  If  the  in¬ 
put  vectors  ■w[n]  are  independent  and  persistent,  and  w[n]  and  Pi[n]  are  mutually 
independent  Gaussian  random  variables  independent  of  all  earlier  values  of  Pi[^])  the 
LMS  algorithm  can  be  shown  to  converge  in  the  mean  to  q^o,  the  Wiener  estimate 
of  qi[n]  as  n  ^  CO.  Even  when  these  stringent  conditions  are  not  satisfied,  the  LMS 
algorithm  is  often  still  effective  and  used  in  practice  because  of  its  computational 
simplicity  compared  to  algorithms  such  as  those  described  in  the  previous  two  sec- 
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tions.  One  problem  with  the  LMS  algorithm  given  by  Eq.  2.40  is  the  norm  of  the 
correction  ti{pi[n]  -  qi[n]w[n])w®[n]  depends  on  the  norm  of  the  input  vector  w[n]. 
This  can  result  in  a  convergence  rate  that  depends  on  the  input  sequence.  To  reduce 
this  dependence,  the  normalized  LMS  [54]  algorithm  uses  ij,  =  p,/{c+  ||w[n]|p).  The 
factor  of  ||w[n]|p  normalizes  the  step  size  by  the  norm  of  the  input  vector,  reduc¬ 
ing  the  dependence  of  the  convergence  rate  on  the  input  sequence.  The  constant  c 
prevents  fj,  from  growing  too  large  if  ||w[n]|p  grows  very  small. 

The  normalized  LMS  algorithm  can  be  used  to  estimate  all  the  rows  of  Q[n] 
simultaneously  by 

Q[n  +  1]  =  Q[n]  -h  (pN  "  QNw[n])  w«[n].  (2.41) 

This  equation  reveals  the  normalized  LMS  estimator  has  the  same  fixed  point  as  the 
CNLLS  estimator.  When  the  estimator  errors  for  the  rows  of  Q[n]  are  orthogonal  to 
the  source  array  vector  w[n]  the  LMS  estimator  will  have  same  convergence  problems 
caused  by  an  impersistent  input  sequence  as  the  CNLLS  estimator.  Comparing 
Eqs.  2.41  and  2.19-2.22  indicates  the  LMS  estimator  is  equivalent  to  the  CNLLS 
estimator  with  #[n]  held  fixed  at  /j1.  Thus,  the  LMS  estimator  is  not  using  any 
information  about  the  past  values  of  w[n]  to  update  Q[n].  The  CNLLS  estimator 
will  generally  converge  and  respond  to  transients  better  because  of  its  use  of  this 
information.  If  the  weight  sequence  w[n]  is  distributed  such  that  the  E  =  I, 
the  CNLLS  and  LMS  estimators  will  be  equivalent  in  the  mean  sense.  However, 
the  constraint  of  exciting  a  fixed  pressure  field  will  generally  require  w[n]  to  remain 
fixed  or  nearly  fixed  once  the  estimate  of  Q[n]  has  converged  close  to  the  true  Q, 
and  thus  generally  E  ^1.  In  this  situation,  the  LMS  estimator  is  not  expected 
to  perform  as  well  as  the  CNLLS  estimator.  The  experiments  presented  in  Chapter 
4  verify  this. 
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Chapter  3 


Simulation  Results 

This  chapter  presents  the  results  of  using  the  control  algorithm  presented  in  the 
last  chapter  in  several  simulated  ocean  environments.  All  of  the  environments  were 
simulated  with  the  finite-element  parabolic  equation  (FEPE)  model  for  acoustic 
propagation  [21].  This  is  an  especially  good  model  to  use  for  these  simulations 
because  it  is  not  a  mode-based  propagation  model.  Consequently,  this  eliminates  any 
concerns  about  the  propagation  model  generating  a  single  mode  because  it  assumes 
the  propagation  takes  place  in  discrete  modes.  The  shallow  water  environments  are 
based  on  observed  sound  speed  profiles  from  the  continental  shelf  in  the  region  of  41° 
N  71°  W  [58].  Figure  3-1  indicates  the  specific  transect  where  the  profiles  used  in  the 
experiments  were  measured.  Section  3.1  confirms  the  algorithm  works  for  a  simple 
range-  and  time-invariant  environment.  Sections  3.2  and  3.3  simulate  two  common 
shallow  water  bathymetric  features:  a  rock  outcropping  and  a  downsloping  wedge. 
Both  of  these  features  affect  acoustic  propagation  and  would  cause  problems  for 
an  open  loop  control  algorithm  which  assumed  the  modes  propagated  adiabatically 
downrange.  The  final  section,  Section  3.4,  simulates  the  propagation  of  a  solitary 
internal  wave  through  the  feedback  volume  over  the  course  of  roughly  twenty-five 
minutes.  The  successful  performance  of  the  feedback  control  algorithm  in  all  these 
environments  offers  encouragement  for  its  successful  deployment  in  an  oceanographic 
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experiment. 


Figure  3-1:  Hydrographic  Transect  Providing  Environmental  Data  for  Simulations 


3.1  Range-Invariant  Environment 

The  simulations  in  this  section  verify  that  the  feedback  control  algorithm  works  for  a 
simple  but  realistic  shallow  water  environment.  In  addition,  the  performance  of  the 
algorithm  is  evaluated  over  a  range  of  estimator  parameters.  For  the  CNLLS  esti¬ 
mator,  the  parameter  was  the  exponential  forgetting  factor,  while  the  Kalman  filter 
trials  varied  the  autoregressive  parameter  a  and  the  innovation  covariance  Pf  used 
by  the  Kalman  filter’s  state  space  model  for  the  replica  matrix  dynamics  in  Eqs.  2.34- 
2.37.  The  environment  used  was  time-invariant  and  homogeneous  in  range  with  the 
vertical  profile  shown  in  Figure  3-2.  The  experiments  were  all  configured  with  the 
source  and  feedback  arrays  separated  by  1  km,  and  the  propagation  frequency  was 
400  Hz.  At  this  frequency,  the  environment  supports  nine  trapped  modes.  All  of  the 
evanescent  modes  and  continuum  energy  will  be  completely  attenuated  over  the  km 
distance  between  the  source  and  feedback  arrays  [18].  The  source  array  consists  of 
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10  omni-directional  point  sources,  and  the  feedback  array  of  19  hydrophones.  The 
replica  matrix  at  this  frequency  for  the  environment  was  computed  using  a  version 
of  FEPE  modified  to  give  complex  replica  vectors  rather  than  transmission  loss. 


Figure  3-2:  Vertical  Profile  of  Range-Invariant  Environment 

Figure  3-3  plots  the  mode  shapes  for  the  nine  trapped  acoustic  modes  of  this 
environment  as  computed  by  FEMODE,  a  finite-difference  mode  computation  pro¬ 
gram.  Figure  3-3  also  displays  the  shape  of  the  sound  speed  profile  at  the  left  edge 
of  the  plot  for  purpose  of  comparison  with  the  mode  shapes.  The  feedback  array 
samples  the  mode  shapes  well,  and  so  has  a  condition  number  of  roughly  1.25. 
Because  is  well-conditioned,  choosing  the  source  array  weights  to  minimize  the 
error  between  the  desired  and  obtained  profiles  in  the  pressure  domain  will  give  a 
high  fidelity  single  mode. 

The  experiments  each  consisted  of  100  independent  trials  attempting  to  excite 
mode  two.  Each  trial  started  from  a  random  Q[0]  and  then  allowed  the  algorithm 
to  converge  for  250  iterations.  The  observation  noise  at  the  feedback  array  was 
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Figure  3-3:  Mode  Shapes  for  Range-Invariant  Environment 

spatially  white,  and  20  dB  below  the  energy  in  the  desired  pressure  profile  for  mode 
two.  This  SNR  is  relatively  high  but  realistic  given  the  short  distance  between  the 
source  and  feedback  arrays.  For  the  scenario  shown  in  Fig.  1-1,  the  excited  field  must 
be  strong  at  the  feedback  array  if  it  is  to  be  measurable  at  the  observation  array, 
which  is  much  further  downrange.  If  the  pressure  field  at  the  feedback  array  is  not 
significantly  above  the  observation  noise,  the  observation  array  may  have  difficulty 
separating  the  signal  due  to  the  mode  input  to  the  observation  volume  from  noise. 

Figure  3-4  shows  the  convergence  of  the  control  algorithm  using  the  CNLLS 
estimator  for  Q.  The  figure  plots  the  mean  performance  over  100  trials  for  three 
different  forgetting  factors,  as  well  as  the  least  squares  bound  obtainable  from  perfect 
knowledge  of  Q.  The  three  forgetting  factors  correspond  to  effective  window  lengths 
of  50,  100  and  200  iterations,  where  the  effective  window  length  Nj  for  a  forgetting 
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factor  7  is  defined  to  be  the  number  of  iterations  required  for  it  to  decay  to  0.1,  i.e., 

=  0.1. 


The  curves  in  Figure  3-4  plot  the  mean  signal-to-error  ratios  (SERs)  averaged  over 
all  trials  for  each  experiment.  The  SER  is  defined  to  be  the  ratio  of  the  energy  in 
the  desired  pressure  profile  to  the  energy  in  the  error  at  each  iteration,  i.e.. 


10  logio 
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where  p[n]  is  the  pressure  observed  at  the  feedback  array  for  the  current  iteration 
including  all  observation  noise.  The  error  term  ||prf-p[n]|p  appears  in  the  denomina¬ 
tor  of  the  SER  expression,  while  the  numerator  UprflP  is  fixed  for  these  experiments. 
Thus,  decreasing  the  error  between  the  desired  and  observed  pressure  vectors  will  in¬ 
crease  the  SER.  The  experiments  indicate  that  the  performance  does  not  vary  much 
with  changes  in  forgetting  factor.  The  “LS  Bound”  at  55  dB  is  the  least-squares 
bound  on  the  performance  that  could  be  attained  with  perfect  knowledge  of  Q  and 
no  observation  noise.  This  bound  is  finite  because  of  numerical  differences  between 
the  methods  for  computing  the  pressure  profiles  for  mode  two  and  the  replica  matrix 
Q.  As  a  result,  a  small  component  of  p^  is  not  in  the  span  of  Q. 

Because  the  observed  pressure  includes  noise,  there  is  a  practical  bound  on  the 
curves  in  Fig.  3-4  considerably  below  the  LS  bound.  Even  if  the  algorithm  perfectly 
excited  the  desired  pressure  field,  the  observation  noise  would  still  limit  the  per¬ 
formance  to  20  dB.  Figure  3-5  shows  the  results  of  several  experiments  using  the 
CNLLS  estimator  with  decreasing  observation  noise  (increasing  SNR).  Again,  each 
curve  is  the  average  performance  over  100  trials  of  the  same  experiment.  In  each 
case,  the  SER  approaches  the  SNR  demonstrating  the  performance  of  the  algorithm 
in  this  environment  is  limited  by  the  level  of  the  observation  noise.  Figure  3-6  shows 
the  desired  pressure  profile  at  the  feedback  array  against  the  pressure  obtained  by  a 
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typical  trial  of  this  experiment  with  effective  window  length  100.  While  the  pressure 
obtained  is  not  a  perfect  fit  to  the  desired  profile,  it  is  a  good  approximation. 


Figure  3-4:  Ratio  of  Energy  in  Desired  Profile  to  Energy  in  Mean  Error  for  CNLLS 
Estimator 

Figures  3-7  through  3-9  show  the  results  of  several  simulations  using  the  Kalman 
filter  channel  estimator  with  different  values  for  a  and  Pf  in  the  Kalman  filter’s  state 
space  model  for  the  dynamics  of  Q.  Note  that  this  is  only  changing  the  parameters 
of  the  model  used  by  the  estimator,  and  does  not  change  the  properties  of  the  real 
channel,  which  is  time  invariant.  These  experiments  were  performed  with  the  same 
configuration,  initialization  and  number  of  trials  as  the  CNLLS  experiments.  For 
a  =  0.98  and  0.99,  the  experiments  indicate  the  control  algorithm  performed  best 
when  a  relatively  large  innovation  covariance  Pf  was  used  in  the  Gauss-Markov 
model  for  the  evolution  of  the  rows  of  Q  i.e.,  Eq.  2.26.  Choosing  Pf  =  10“®I  gave 
the  best  performance  for  these  values  of  a.  The  mean  absolute  value  of  the  elements 
of  the  true  Q  is  1.2  x  10“^,  so  the  standard  deviation  of  the  innovation  process  for 
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igure  3-5:  SER  for  CNLLS  Estimator  at  different  SNR  levels 


Figure  3-6:  Vertical  Profile  of  Pressure  for  Typical  Trial  of  CNLLS  Estimator 


these  elements  is  on  the  same  order  as  the  mean  replica  matrix  element.  When 
a  =  1.0,  Figure  3-9  shows  the  estimator  with  Pf  =  10~®I  does  noticeably  worse 
than  the  others.  Solving  for  the  steady-state  error  covariance  of  the  Kalman  filter 
reveals  that  for  a  =  1.0  the  error  covariance  of  the  Kalman  filter  is  much  larger  for 
Pf  =  10“®I.  It  is  sensible  that  a  large  value  of  Pf  would  make  the  error  variance  grow 
even  faster  and  cause  worse  performance  than  smaller  values  of  Pf  when  a  =  1.0  so 
there  is  no  decay.  For  all  three  values  of  a  shown  here,  experiments  using  even  larger 
values  for  Pf  caused  worse  performance.  This  is  reasonable  because  a  larger  Pf  would 
mean  the  innovation  portion  of  the  Gauss-Markov  process  f  would  overwhelm  the 
true  values  of  the  time-invariant  Q  the  algorithm  is  attempting  to  estimate.  Figure 
3-10  plots  the  desired  and  obtained  vertical  pressure  profiles  for  a  typical  trial  with 
a  =  0.99  and  Pf  =  10~®I.  Similar  to  the  result  shown  in  Figure  3-6,  the  profile 
obtained  by  the  control  algorithm  is  a  close  but  not  perfect  fit  to  the  desired  profile. 

Figure  3-11  plots  the  results  of  a  simulation  comparing  the  CNLLS  with  7=1 
to  the  Kalman  filter  with  a  =  1  and  Pf  =  0.  It  can  be  demonstrated  that  for 
these  parameter  choices,  the  estimators  are  equivalent  if  the  condition  number  is 
not  limited  in  the  CNLLS  estimator.  The  curves  in  Fig.  3-11  are  very  close  but 
not  identical.  Two  factors  may  have  caused  this  discrepancy.  First,  each  curve  is 
an  average  of  100  trials  with  the  appropriate  estimator.  Each  trial  started  with  a 
different  randomly  chosen  Q[0].  In  the  limit  of  an  infinite  set  of  trials,  these  averages 
should  agree,  but  the  sample  averages  for  each  estimator  need  not  be  identical  after 
100  trials.  Second,  the  CNLLS  estimator  did  limit  the  condition  number  of  #“^[n], 
so  the  estimators  are  not  exactly  equivalent  when  $  becomes  poorly  conditioned. 
The  increasing  difference  between  the  curves  with  increasing  n  supports  this  expla¬ 
nation.  If  the  condition  number-limiting  is  a  factor,  it  will  only  matter  once  w[n]  has 
converged  and  #[n]  becomes  poorly-conditioned.  This  will  take  several  iterations, 
and  the  effect  will  become  more  severe  as  $  becomes  more  poorly  conditioned  over 
time.  The  curves  in  Fig.  3-11  start  very  close  together  and  only  diverge  slowly  as  the 
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would  be  expected  as  the  conditioning  of  ^  grew  worse.  This  supports  the  hypothe¬ 
sis  the  condition  number-limiting  is  the  primary  factor  responsible  for  the  difference 
between  the  estimators.  On  the  whole,  the  close  agreement  between  the  estimators 
in  this  experiment  offers  reassurance  that  the  estimators  are  correctly  implemented. 


Figure  3-7;  Ratio  of  Energy  in  Desired  Profile  to  Energy  in  Mean  Error  for  Kalman 
Filter  Estimator  with  a  =  0.98 


3.2  Rock  Outcropping 

This  section  examines  the  performance  of  the  control  algorithm  when  the  region  be¬ 
tween  the  source  and  feedback  arrays  includes  a  rock  outcropping.  Such  features  are 
common  in  shallow  water  environments,  and  the  strong  interaction  of  the  pressure 
field  with  the  bottom  means  these  bathymetric  features  have  an  important  effect 
on  propagation.  Given  very  accurate  and  detailed  environmental  measurements,  an 
open  loop  control  algorithm  could  theoretically  still  generate  a  single  mode  in  the 
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Figure  3-8:  Ratio  of  Energy  in  Desired  Profile  to  Energy  in  Mean  Error  for  Kalman 
Filter  Estimator  with  a  =  0.99 


Figure  3-9:  Ratio  of  Energy  in  Desired  Profile  to  Energy  in  Mean  Error  for  Kalman 
Filter  Estimator  with  a  =  1.0 


Figure  3-10:  Vertical  Profile  of  Pressure  for  Typical  Trial  of  Kalman  Filter  Estimator 
with  a  =  0.99,  Pf  =  10“®I. 


Figure  3-11:  Comparison  of  CNLLS  and  Kalman  filter  estimators 
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presence  of  such  an  outcropping.  An  echosounder  survey  could  give  a  reasonably  ac¬ 
curate  bathymetric  profile  of  the  environment,  though  it  would  still  have  uncertainty 
in  depth  on  the  order  of  a  quarter  to  half  wavelength  (l-2m).  The  acoustic  effect 
of  this  uncertainty  would  probably  be  small,  but  not  insignificant.  Measurements 
of  the  geophysical  parameters  for  the  specific  location  would  require  exhaustive  core 
and  sediment  sampling.  However,  representative  measurements  of  these  parameters 
exist  for  most  coast  regions,  and  generally  these  measurements  would  be  acceptable 
approximations.  The  open  loop  algorithm  would  start  with  the  available  environ¬ 
mental  data  and  try  to  compute  the  replica  matrix  from  acoustic  propagation  models. 
The  source  array  weights  would  then  be  set  to  be  those  the  computed  replica  ma¬ 
trix  predicted  would  give  the  desired  pressure  field  at  the  start  of  the  observation 
volume.  In  contrast,  the  feedback  control  algorithm  estimates  the  replica  matrix  di¬ 
rectly  from  acoustic  observations,  removing  the  intermediate  step  of  computing  the 
replica  matrix  from  hydrographic  and  bathymetric  observations.  To  whatever  extent 
the  geophysical  parameters  of  the  environment  are  relevant  to  the  mode  excitation 
experiment,  they  should  manifest  themselves  in  the  observed  pressure  data. 

Figure  3-12  shows  the  propagation  environment  with  the  outcropping.  Other 
than  the  6  m  tall  outcropping  at  r  =  200  m,  this  environment  is  identical  to  the 
range-invariant  environment  used  in  Section  3.1.  The  source  array  is  at  r  =  0  m  and 
the  feedback  array  at  r  =  1000  m.  The  outcropping  has  physical  properties  similar 
to  those  of  limestone  or  basalt  [32].  The  experiments  in  this  section  will  assume 
the  extreme  scenario  that  no  a  priori  information  was  available  about  the  rock  out¬ 
cropping.  The  open  loop  algorithm  assumes  it  is  operating  in  the  range-invariant 
environment  shown  in  Fig.  3-2.  Generally,  better  information  will  be  available,  but 
this  exaggerated  mismatch  demonstrates  the  attractive  possibility  of  using  the  feed¬ 
back  array  without  any  a  priori  knowledge  of  the  environment. 

The  presence  of  this  outcropping  has  a  strong  effect  on  the  acoustic  propagation. 
Figure  3-13  plots  the  performance  of  the  feedback  control  algorithm  using  both 
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Figure  3-12:  Rock  Outcropping  Environment 

estimators,  as  well  as  for  an  open-loop  controller  which  assumes  the  environment  is 
homogeneous  in  range  with  the  vertical  profile  observed  at  the  source  array.  All  the 
curves  are  the  average  over  100  trials  initialized  with  random  replica  matrices.  The 
CNLLS  estimator  used  an  effective  window  length  of  100,  while  the  Kalman  filter 
had  a  =  0.99  and  Pf  =  10~®I.  The  figure  shows  both  estimators  converge  to  allow 
the  feedback  controller  to  excite  a  single  mode.  However,  the  very  poor  performance 
of  the  open  loop  controller  indicates  the  importance  of  accounting  for  the  effect  of 
the  outcropping  on  propagation.  By  assuming  the  environment  was  range-invariant 
and  that  modes  would  propagate  adiabatically,  the  open-loop  controller  generated  a 
pressure  profile  in  the  far  field  which  was  a  poor  approximation  to  the  desired  pressure 
profile,  and  likely  unacceptable  as  a  source  signal  for  oceanographic  observations. 
Figure  3-14  shows  vertical  pressure  profiles  at  the  feedback  array  for  one  trial  after 
the  estimators  had  converged.  From  these  plots,  it  is  apparent  that  the  feedback 
control  algorithm  gives  a  good  approximation  to  the  desired  pressure  profile  using 
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either  estimator.  The  open-loop  controller  generates  a  profile  that  is  unacceptable 
as  a  single  mode  source. 


Figure  3-13:  Performance  of  Estimators  For  Rock  Outcropping  Environment 

Figures  3-15  and  3-16  present  a  revealing  contrast  between  the  pressure  fields 
excited  by  the  open  loop  and  feedback  control  algorithms.  Both  figures  plot  acoustic 
intensity  as  a  function  of  range  and  depth,  where  high  intensity  is  indicated  by  a 
light  shading.  The  intensity  has  been  normalized  to  remove  the  effect  of  cylindrical 
spreading.  The  desired  mode  was  mode  two  and  the  environment  was  the  rock 
outcropping  environment  shown  in  Figure  3-12  in  both  cases.  Figure  3-15  shows  the 
pressure  profile  versus  range  and  depth  for  the  field  excited  by  the  open  loop  control 
algorithm,  as  well  as  the  vertical  profile  of  the  magnitude  of  the  desired  pressure  field 
for  mode  two  on  the  left  side  of  the  figure  for  reference.  Close  to  the  source  array, 
the  pressure  field  can  be  seen  to  be  settling  towards  mode  two,  with  nulls  developing 
at  the  surface  and  2:  =  24  m.  When  the  wave  is  incident  upon  the  outcropping  at 
r  =  200  m,  the  energy  is  scattered  into  several  modes,  generating  the  complicated 
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Figure  3-14:  Vertical  Profile  of  Pressure  for  a  Typical  Trial  in  Rock  Outcropping 
Environment 

interference  pattern  seen  downrange  of  the  outcropping.  Based  on  this  pattern,  it 
is  clear  the  pressure  field  entering  the  observation  volume  consists  of  several  modes, 
and  not  the  desired  single  mode. 

Figure  3-16  plots  the  intensity  for  the  feedback  controller  using  the  CNLLS  es¬ 
timator.  For  this  algorithm,  the  source  array  initially  excites  a  very  complicated 
interference  pattern  containing  a  rich  variety  of  modes.  When  this  pressure  field 
impinges  upon  the  rock  outcropping,  the  energy  scatters  such  that  all  the  undesired 
modes  cancel  out,  leaving  only  mode  two  propagating.  By  the  time  the  field  reaches 
the  feedback  array  at  the  start  of  the  observation  volume  (r  =  1000  m),  there  is  a 
well-developed  pressure  null  at  a  depth  of  2:  =  24  m.  This  corresponds  well  with  the 
desired  vertical  pressure  profile  shown  at  the  left  side  of  the  figure.  The  feedback 
control  algorithm  successfully  estimates  and  compensates  for  the  scattering  effect  of 
the  rock  outcropping  to  produce  the  desired  pressure  field  in  the  far  field. 
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These  results  can  be  interpreted  in  terms  of  the  matched  signal  paradigm  dis¬ 
cussed  in  Section  1.5.  The  open  loop  algorithm  assumes  an  adiabatic  environment 
and  transmits  the  signal  desired  at  the  feedback  array  immediately  from  the  source 
array.  The  outcropping  distorts  the  transmitted  mode  two,  coupling  the  energy  into 
several  modes.  The  feedback  control  algorithm  estimates  the  environmental  effect 
on  the  propagation,  then  transmits  a  complicated  waveform  which  is  matched  to 
the  environment  so  the  pressure  field  at  the  start  of  the  observation  volume  consists 
solely  of  mode  two.  This  is  analogous  to  the  work  done  by  Parvulescu  and  Clay  [46], 
[45]  only  here  the  desired  waveform  is  concentrated  in  the  mode  domain  rather  the 
time  domain.  The  feedback  control  algorithm  transmits  a  waveform  that  is  initially 
not  concentrated  in  the  mode  domain  such  that  the  pressure  field  at  the  start  of 
the  observation  volume  contains  a  single  mode.  The  open  loop  controller  does  not 
account  for  any  channel  effects,  and  the  field  it  generates  at  the  start  of  the  obser¬ 
vation  volume  is  not  concentrated  in  the  mode  domain,  but  has  energy  in  several 
modes. 

3.3  Downsloping  Wedge 

A  downsloping  wedge  is  a  common  feature  in  many  shallow  water  environments. 
This  section  presents  the  results  of  simulating  the  deployment  of  the  feedback  con¬ 
trol  algorithm  in  a  2°  downsloping  wedge  using  the  profiles  observed  in  [58].  For  the 
portions  of  the  wedge  deeper  than  the  bottom  used  in  the  previous  section,  the  deep¬ 
est  measurement  of  sound  speed  from  the  previous  section  is  extended  downward  to 
form  a  isovelocity  layer  above  the  bottom  in  the  wedge.  The  experiments  described 
in  this  section  will  compare  the  performance  of  the  feedback  control  algorithm  to 
the  open-loop  algorithm,  which  assumes  the  modes  excited  at  the  source  array  prop¬ 
agate  downslope  without  coupling.  In  reality,  the  slope  causes  coupling  among  the 
modes  as  they  propagate  downrange  so  the  field  excited  by  the  open-loop  algorithm 
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contains  several  undesired  modes  when  it  enters  the  observation  volume.  The  feed¬ 
back  control  algorithm  successfully  compensates  for  the  coupling  introduced  by  the 
sloping  bottom  to  excite  the  desired  pressure  field  at  the  start  of  the  observation 
volume. 


Figure  3-17;  Downsloping  Wedge  Environment. 


Figure  3-17  shows  the  environment  used  for  this  simulation.  The  ten  element 
source  array  is  located  at  r  =  0,  and  the  nineteen  hydrophone  feedback  array  is  at 
r  =  1000  m.  The  depths  at  the  source  and  feedback  array  locations  are  30.75  m 
and  65.75  m,  respectively.  For  the  transmission  frequency  of  400  Hz,  the  channel 
supports  eight  modes  at  the  source  array  and  nineteen  modes  at  the  feedback  array. 
The  depth  at  the  source  array  is  only  slightly  shallower  than  the  cutoff  depth  for 
mode  nine.  Mode  nine  is  initially  an  evanescent  mode,  but  as  the  channel  becomes 
deeper,  it  becomes  trapped  before  the  energy  propagating  at  that  wavenumber  has 
been  significantly  attenuated.  This  phenomenon,  known  as  mode  capture,  has  been 
observed  both  in  laboratory  experiments  [31]  and  numerical  simulations  ,  [60].  The 
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simulations  in  this  section  verify  that  the  replica  matrix  model  for  propagation  used 
by  the  feedback  control  algorithm  incorporates  continuum  effects  in  propagation 
sufficiently  to  model  mode  capture  in  downslope  propagation.  The  sound  speed 
profile  and  other  environmental  parameters  are  the  same  as  those  used  in  Section  3.2. 

One  hundred  trials  of  the  feedback  control  algorithm  were  run  using  both  the 
Kalman  filter  and  the  CNLLS  estimators.  The  Kalman  filter  estimator  used  Pf  = 
10“®I  and  a  =  0.99  and  the  CNLLS  used  a  forgetting  factor  giving  an  effective 
window  length  of  100  iterations.  Each  trial  was  initialized  with  a  random  Q[0]  and 
was  allowed  to  run  for  250  iterations  attempting  to  excite  mode  one.  The  observation 
noise  was  spatially  white  and  20  dB  below  the  energy  of  the  desired  signal,  p^.  Figure 
3-17  shows  the  desired  pressure  field  shape  as  a  dotted  line  near  the  feedback  array. 
Figure  3-18  shows  the  mean  SER  for  these  trials,  as  well  as  the  performance  of  the 
open  loop  algorithm.  The  open  loop  algorithm  chooses  the  source  array  weights  to 
excite  mode  1  assuming  the  mode  would  propagate  adiabatically  downslope  from  the 
source  array.  The  figure  shows  that  both  the  Kalman  filter  and  CNLLS  estimators 
work  well,  but  the  error  in  the  open  loop  controller  exceeds  the  energy  of  the  desired 
signal.  Figure  3-19  plots  the  vertical  profile  of  pressure  magnitude  for  typical  trials 
using  each  estimator  and  the  open  loop  algorithm.  Regardless  of  which  channel 
estimator  is  used,  the  pressure  field  generated  by  the  feedback  controller  is  a  closer 
approximation  to  the  desired  pressure  profile  than  the  field  generated  by  the  open 
loop  controller. 

3.4  Solitary  Internal  Wave 

This  section  examines  the  performance  of  the  feedback  control  algorithm  in  a  simu¬ 
lated  environment  with  a  propagating  solitary  internal  wave.  The  previous  simula¬ 
tions  in  this  chapter  studied  environments  where  time-invariant  mode  coupling  was 
induced  by  bathymetric  features  such  as  a  sloping  bottom  or  rock  outcropping.  One 
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Figure  3-18:  Mean  SER  Convergence  for  Downsloping  Wedge  Environment 


Figure  3-19:  Typical  Vertical  Pressure  Profiles  for  Downsloping  Wedge  Environment 


of  the  primary  motivations  for  developing  a  single  mode  source  such  as  the  feedback 
control  algorithm  is  to  study  the  acoustic  coupling  introduced  by  internal  waves  in 
the  observation  volume.  If  the  feedback  control  algorithm  is  to  function  successfully 
as  a  single  mode  source  in  such  an  environment,  it  must  be  capable  of  exciting  a 
single  mode  even  when  solitary  waves  are  propagating  through  the  feedback  volume. 
This  simulation  models  the  time- varying  acoustic  mode  coupling  between  the  source 
and  feedback  arrays  introduced  by  a  solitary  internal  wave  propagating  through  the 
feedback  volume.  Based  on  studies  of  internal  waves  supported  on  density  gradients 
in  coastal  regions  [61]  [62]  [63],  the  thermocline  at  2  =  10  m  of  the  simulated  envi¬ 
ronment  is  displaced  by  10sech^((r-  1300  +  t)A:i„t)  m  for  1000  <r  <  1600,  where  the 
horizontal  wavenumber  of  the  internal  wave  is  ki^t  —  (27r/300)m“^  Starting  with  a 
typical  sound  speed  profile  observed  in  [58] ,  this  displacement  profile  gives  the  sound 
speed  profile  shown  in  Fig.  3-20.  In  this  figure,  lighter  regions  indicate  water  with  a 
higher  sound  speed,  while  darker  regions  indicate  colder  water  with  a  slower  sound 
speed. 

The  source  array  is  located  at  r  =  0  m,  and  the  feedback  array  at  r  =  1000  m.  The 
solitary  wave  propagates  towards  r  =  0  m  with  a  velocity  of  1  m/s,  consistent  with 
the  velocities  observed  for  these  waves  in  [61].  The  replica  matrices  were  computed 
every  10  m  (or  equivalently  every  10  seconds)  using  FEPE.  The  frequency  is  400  Hz, 
which  gives  nine  trapped  modes  for  the  water  depth  of  34  m.  The  simulation  used  the 
same  bottom  parameters  and  source  and  feedback  arrays  as  Section  3.1.  The  feedback 
control  algorithm  attempts  to  excite  mode  two  starting  from  a  randomly  chosen 
Q[0]  and  using  an  iteration  time  of  1  sec.  The  slowest  mode  in  this  environment 
propagates  from  the  source  array  to  the  feedback  array  in  roughly  0.8  sec,  so  this 
choice  of  iteration  time  insures  the  field  observed  at  the  feedback  array  reflects  the 
most  recent  source  array  weights.  The  observation  noise  is  again  spatially  white  and 
20  dB  below  the  pressure  profile  for  mode  two. 

The  simulation  was  run  with  both  the  CNLLS  and  Kalman  Filter  channel  estima- 
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Figure  3-20:  Sound  speed  profile  c{r,  z)  for  solitary  internal  wave  with  sech^  profile 


tors.  The  former  used  a  forgetting  factor  of  7  =  0.9772,  giving  an  effective  window 
length  of  100  iterations.  The  Kalman  filter  estimator  set  a  =  0.99  and  Pf  =  10“®I. 
Figures  3-21  and  3-22  show  the  performance  of  the  feedback  control  algorithm  using 
the  CNLLS  and  Kalman  filter  estimators,  respectively.  Because  the  internal  wave 
is  propagating  with  velocity  1  m/s,  the  time  axes  of  Figs.  3-21  and  3-22  can  also 
be  interpreted  as  the  initial  range  location  in  Fig.  3-20  of  the  sound  speed  profile 
currently  located  at  the  source  array.  Each  figure  also  shows  the  performance  of  the 
open  loop  algorithm  in  the  same  environment. 

Figures  3-21  and  3-22  clearly  demonstrate  the  feedback  control  algorithm  using 
either  estimator  performs  better  than  the  open  loop  algorithm  in  exciting  the  desired 
mode.  The  open  loop  control  algorithm  observes  the  initial  sound  speed  profile  at 
the  source  array,  and  assumes  the  feedback  volume  is  range-invariant  with  this  profile 
throughout  the  experiment.  Initially,  this  assumption  is  valid  and  attains  the  limit 
on  SER  imposed  by  the  observation  noise.  Once  the  solitary  wave  starts  entering 
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Figure  3-21:  Performance  of  feedback  control  algorithm  using  CNLLS  estimator 
compared  with  open  loop  controller  for  solitary  wave  propagating  through  feedback 
volume. 


Time  (sec) 

Figure  3-22:  Performance  of  feedback  control  using  Kalman  filter  estimator  compared 
with  open  loop  controller  for  solitary  wave  propagating  through  the  feedback  volume. 
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the  feedback  volume  the  open  loop  algorithm’s  performance  deteriorates  severely  and 
does  not  recover  until  the  wave  has  largely  propagated  past  the  source  array.  The 
Kalman  filter,  in  particular,  tracks  the  time-varying  coupling  introduced  by  the  soli¬ 
tary  wave  to  maintain  a  consistently  high  fidelity  mode  after  an  initial  transient.  The 
pressure  at  the  feedback  array  for  the  CNLLS  has  significant  oscillations,  making  its 
pressure  field  a  less  desirable  source  for  downrange  observers  than  the  field  generated 
by  the  Kalman  filter.  The  CNLLS  algorithm  assumes  the  dynamics  of  the  rows  of 
the  replica  matrix  are  stationary  over  the  length  of  the  eflFective  window  introduced 
by  the  forgetting  factor.  The  deterministic  propagation  of  the  solitary  wave  through 
the  feedback  volume  violates  this  assumption.  Simulations  run  with  time-invariant 
environments  chosen  from  the  replica  matrices  representing  the  solitary  wave  inside 
the  feedback  volume  show  the  CNLLS  estimator  converges  under  these  conditions. 
This  also  suggests  the  oscillations  are  caused  by  the  nonstationary  environment.  The 
Kalman  filter  does  not  assume  the  replica  matrix  is  stationary,  and  is  able  to  ex¬ 
cite  a  high  fidelity  single  mode  while  a  strong  solitary  wave  propagates  through  the 
feedback  volume  inducing  time- varying  coupling. 

3.5  Summary 

This  chapter  has  demonstrated  that  the  feedback  control  algorithm  is  a  viable 
method  for  exciting  a  single  mode  in  several  simulated  coastal  environments  ranging 
from  a  simple  adiabatic  range-  and  time-invariant  channel  to  a  complicated  range- 
and  time- varying  propagation  environment  generated  by  a  propagating  solitary  inter¬ 
nal  wave.  These  encouraging  results  suggest  that  the  next  logical  step  is  to  evaluate 
the  performance  of  the  feedback  control  algorithm  in  a  scale  model  experiment.  Such 
an  experiment  will  introduce  the  limitations  of  real-time  hardware  and  the  vagaries 
of  physical  acoustic  propagation,  testing  the  robustness  of  the  algorithm.  The  next 
chapter  describes  a  set  of  experiments  which  study  the  feedback  control  algorithm 
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under  these  conditions. 
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Chapter  4 


Laboratory  Waveguide 
Experiments 


This  chapter  describes  a  series  of  experiments  using  the  feedback  control  algorithm 
to  control  the  pressure  field  in  a  laboratory  waveguide.  These  experiments  investi¬ 
gated  the  performance  of  the  feedback  control  algorithm  implemented  on  real-time 
hardware  and  with  actual  rather  than  simulated  acoustic  propagation.  Section  4.1 
describes  the  laboratory  waveguide  used  for  these  experiments  as  well  as  the  trans¬ 
ducers  and  signal  processing  hardware.  The  acoustic  characteristics  of  the  waveguide 
are  discussed  in  Section  4.2.  The  first  set  of  experiments  confirmed  that  the  feed¬ 
back  control  algorithm  could  successfully  control  the  pressure  field  in  the  waveguide 
without  a  priori  knowledge  of  the  propagation  characteristics.  Section  4.3  describes 
these  experiments  as  well  as  a  second  series  of  experiments  that  explored  which 
parameters  of  the  channel  estimators  gave  the  best  performance  in  the  waveguide. 
Section  4.4  presents  the  results  of  a  series  of  experiments  investigating  the  transient 
response  of  the  feedback  algorithm  when  an  acoustic  scatterer  was  introduced  into 
the  waveguide  after  the  algorithm  had  been  allowed  to  converge. 
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4.1  Waveguide  and  Equipment  Description 


The  laboratory  waveguide,  or  flume,  used  for  these  experiments  was  originally  de¬ 
signed  for  fluid  dynamics  rather  than  acoustics  experiments.  Consequently,  the  flume 
is  not  an  ideal  acoustic  waveguide,  and  little  information  is  available  about  the  as¬ 
pects  of  construction  germane  to  its  acoustic  properties.  Much  of  what  follows  about 
the  materials  or  structure  is  conjecture  based  on  inspection  of  the  flume,  since  no 
records  or  plans  exist  from  the  waveguide’s  construction.  The  tank  is  roughly  20  m 
long,  1  m  deep  and  1.2  m  wide  from  front  to  back.  The  front  face  is  constructed  of 
plates  of  3/8-inch  thick  plexiglass,  joined  by  aluminum  brackets  to  each  other  and 
10  cm  thick  concrete  pilings  roughly  every  2  m.  Fortunately,  the  distance  the  pilings 
project  into  the  flume  is  insigniflcant  compared  to  wavelength  of  these  experiments. 
The  back  wall  appears  to  be  constructed  of  concrete  roughly  10  cm  thick.  The  back 
wall  and  bottom  of  the  flume  were  thickly  painted  or  coated  with  a  latex-type  ma¬ 
terial  in  a  manner  that  left  large  (20-30  cm)  bubbles  irregularly  spaced  in  several 
locations  on  the  bottom.  The  bottom  appeared  to  be  largely  solid  concrete  but  a 
strip  of  signiflcant  width  (30-40  cm)  sounded  hollow,  perhaps  flberglass  over  a  chan¬ 
nel  containing  plumbing  for  the  pumps  used  to  generate  hydrodynamic  flows  in  the 
tank. 

Two  baffling  plates  were  installed  at  either  end  of  the  tank  for  these  experiments 
to  reduce  the  amount  of  energy  reflected,  making  it  a  closer  approximation  to  an 
inflnitely  long  waveguide.  Each  baffler  was  a  steel  plate  covered  with  a  layer  of  VRC 
PB  1-94  anechoic  tiles  manufactured  by  Vector  Research  Corporation  under  con¬ 
tract  from  Seaward  International  and  generously  loaned  to  us  for  these  experiments. 
This  material  consists  of  polyurethane  encapsulated  gaseous  microspheres,  and  was 
molded  into  pyramids  roughly  2  cm  square  at  the  base  and  extending  6  cm  out  from 
the  plate.  The  material  was  designed  for  20  kHz,  rather  than  the  8  kHz  operat¬ 
ing  frequency  used  in  these  experiments.  Preliminary  calibration  data  from  Vector 
Research  indicates  even  at  these  lower  frequencies  the  panels  should  still  provide 
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roughly  10  dB  of  echo  reduction  at  this  lower  frequency.  The  flume  was  fllled  with 
fresh  water  to  a  depth  of  66  cm  at  the  source  array,  which  gave  a  depth  of  69  cm 
at  the  feedback  array  due  to  a  slight  slope  in  the  bottom  of  the  flume.  Assuming 
a  frequency  of  400  Hz  in  the  ocean  experiment,  scaling  these  depths  by  the  same 
factor  of  twenty  as  the  frequency  indicates  this  would  be  equivalent  to  a  13-14  m 
water  depth  for  an  ocean  experiment. 

The  source  array  consisted  of  six  Channel  Industries  3013  ceramic  cylinder  8- 
15  kHz  acoustic  sources,  spaced  10  cm  apart  and  offset  5  cm  from  the  surface  of 
the  water.  The  feedback  array  consisted  of  seven  Benthos  AQ2TS  hydrophones 
spaced  every  8  cm.  The  hydrophones,  along  with  their  pre-amplifiers,  were  sealed 
in  a  polyurethane  tube  fllled  with  mineral  oil  to  minimize  the  acoustic  impedance 
mismatch  with  the  surrounding  water.  A  tapered  lead  weight  ballasted  the  base  of 
the  feedback  array,  allowing  it  to  be  tensioned  to  remain  vertical.  The  pre-amplifiers 
were  calibrated  before  insertion  into  the  array  so  that  all  the  pre-amplifiers  had 
gains  of  36.5  ±  0.1  dB  per  WHOI  standard  [64],  Figure  4-1  provides  a  schematic  of 
the  physical  layout  of  the  waveguide  for  the  experiments.  The  arrays  were  centered 
width-wise  in  the  waveguide,  roughly  60  cm  from  the  front  and  back  walls.  If  the 
source  array  is  defined  as  the  origin  in  range,  with  range  increasing  towards  the 
feedback  array,  the  baffling  plates  were  located  at  r  =  —1.6  m  and  r  =  12.95  m.  The 
feedback  array  was  located  at  r  =  10.9  m. 

The  feedback  control  algorithm  was  implemented  in  the  Acoustic  Modem  Sys¬ 
tem  (AMS-2.0),  a  MATLAB-like  real-time  signal  processing  environment  for  the 
Texas  Instrument  TMS320C40  DSP  chip.  AMS  was  developed  by  Mark  Johnson 
and  Matt  Grund  of  the  Acoustic  Telemetry  Group  at  the  Woods  Hole  Oceanographic 
Institution  to  assist  rapid  development  of  real-time  signal  processing  applications  in 
underwater  acoustics  [65].  The  A/D  and  D/A  conversion  required  for  measuring 
the  pressure  at  the  feedback  array  and  transmitting  the  source  array  weights  were 
handled  by  two  of  the  modems  also  developed  by  the  Acoustic  Telemetry  Group. 
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Figure  4-1:  Configuration  of  Laboratory  Waveguide  Experiments 


One  modem  sampled  the  pressures  observed  at  the  hydrophone  array,  while  the  sec¬ 
ond  modulated  the  complex  source  array  weights  onto  narrowband  sinusoids.  These 
sinusoids  were  fed  to  six  independent  power  amplifiers,  which  drove  the  source  array. 

4.2  Acoustic  Properties  of  the  Laboratory  Waveg¬ 
uide 

This  section  describes  the  acoustic  properties  of  the  laboratory  waveguide  for  the 
experimental  configuration  outlined  in  Section  4.1.  The  water  depth  in  the  waveguide 
was  66  cm  at  the  source  array  and  the  frequency  was  8  kHz.  The  temperature 
of  the  fresh  water  used  to  fill  the  waveguide  was  consistently  25°C,  which  gives  a 
sound  speed  of  1490  m/s.  The  ideal  isovelocity  constant-depth  waveguide  with  a 
pressure-release  surface  and  rigid  bottom  [22]  supports  six  trapped  modes  for  this 
water  depth  and  frequency.  The  bottom  of  the  fiume  is  not  a  homogeneous  perfect 
acoustic  reflector,  making  the  flume  unlikely  to  be  a  good  approximation  to  this 
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ideal  waveguide.  In  addition,  the  boundary  conditions  imposed  across  the  width  of 
the  flume  by  the  plexiglass  front  wall  and  back  wall  of  unknown  composition  further 
complicate  the  acoustic  propagation.  Since  the  flume  has  a  rectilinear  geometry,  the 
width  and  depth  eigenfunctions  of  the  waveguide  should  be  separable.  The  width 
eigenfunctions  should  superimpose  at  the  location  of  the  vertical  feedback  array  to 
form  modes  with  the  vertical  profiles  of  the  depth  eigenfunctions. 

Preliminary  experiments  verified  that  the  acoustic  propagation  between  the  source 
array  and  feedback  array  was  well  modeled  by  an  LTI  system.  Sequential  excitation 
and  scaling  of  source  array  elements  confirmed  the  channel  satisfied  linearity.  Long 
term  experiments  using  the  same  set  of  source  array  weights  indicated  the  channel 
was  effectively  time-invariant,  with  a  correlation  time  on  the  order  of  hours  when 
the  water  in  the  waveguide  was  calm.  These  experiments  also  indicated  that  the  ob¬ 
servation  noise  present  in  the  system  was  more  than  40  dB  below  the  desired  signal 
levels. 

The  sequential  excitation  experiments  also  provided  an  estimate  of  the  replica 
matrix  for  the  waveguide.  Computing  the  SVD  of  this  matrix  revealed  that  only  three 
modes  of  the  system  accounted  for  over  99%  of  the  energy  received  at  the  feedback 
array.  Figures  4-2  and  4-3  plot  the  singular  values  (aq)  and  pressure  modes  (U^) 
for  the  first  four  modes  of  the  estimated  replica  matrix  Q.  The  very  high  SNR  in 
the  waveguide  means  this  should  be  an  accurate  estimate  of  the  true  replica  matrix. 
For  the  remainder  of  the  discussion  in  this  chapter,  the  replica  matrix  estimated  by 
sequentially  exciting  the  sources  will  be  considered  to  be  the  true  Q,  and  the  caret  (•) 
will  be  reserved  for  replica  matrices  estimated  by  the  feedback  control  algorithm.  In 
general,  the  columns  of  Uq  will  not  equal  the  sampled  mode  shapes,  but  they  should 
span  the  same  space  as  the  acoustic  modes.  Projecting  the  theoretical  modes  shapes 
for  an  isovelocity  waveguide  with  a  rigid  bottom  onto  the  observed  Uq  will  indicate 
how  well  these  theoretical  modes  fall  in  the  observed  span  of  the  system.  Performing 
this  projection  reveals  a  significant  portion  of  the  energy  in  the  first  three  theoretical 
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modes  falls  outside  the  span  of  the  first  three  columns  of  Uq.  This  indicates  the 
laboratory  waveguide  is  not  a  good  approximation  to  an  isovelocity  waveguide  with 
a  rigid  bottom.  This  is  not  surprising  given  the  unknown  composition  and  irregular 
construction  of  the  bottom  of  the  flume. 

Repeated  estimations  of  Q  indicate  the  tank  is  well-modeled  as  LTI  over  the 
time-scales  considerably  longer  than  the  propagation  time  between  the  source  and 
feedback  arrays.  Although  the  singular  values  do  vary  slowly  over  time,  they  remain 
roughly  the  same  size  as  those  shown  in  Fig.  4-2.  Based  on  the  values  plotted  in 
Fig.  4-2,  three,  or  at  most  four,  modes  propagate  in  the  laboratory  waveguide. 


The  remainder  of  the  experiments  described  in  this  chapter  use  the  columns  of 
Uq,  or  system  modes,  for  the  desired  pressure  profile.  In  a  well-designed  ocean 
experiment,  the  acoustic  modes  fall  within  the  span  of  the  columns  of  Uq,  allowing 
the  feedback  control  algorithm  to  excite  the  desired  mode.  Demonstrating  the  al¬ 
gorithm  is  capable  of  obtaining  pressure  profiles  in  the  numerical  span  of  Q  in  the 
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Figure  4-3:  Pressure  Magnitude  for  System  Modes  of  Laboratory  Waveguide 

flume  should  give  us  confidence  that  the  algorithm  will  be  successful  in  exciting  the 
desired  acoustic  mode  in  a  well-designed  ocean  experiment. 

4.3  Time-Invariant  Waveguide 

The  experiments  in  this  section  demonstrate  the  feedback  control  algorithm  can  suc¬ 
cessfully  excite  the  desired  pressure  profile  in  a  time-invariant  laboratory  waveguide. 
The  initial  Q[0]  is  the  true  Q  with  noise  added  20  dB  below  the  true  values.  Starting 
with  this  initial  Q[0],  the  feedback  control  algorithm  was  allowed  to  run  for  100  iter¬ 
ations  with  a  0.9  second  iteration  time.  The  desired  pressure  profiles  were  the  system 
modes,  i.e.,  the  columns  of  Uq.  The  experiments  were  run  on  a  different  day  with  a 
slightly  different  Q  than  the  replica  matrix  measured  in  the  last  section.  As  a  result, 
the  desired  pressure  profiles  and  relative  energies  were  similar  but  not  identical  to 
those  shown  in  Figs.  4-3  and  4-2.  The  gains  of  the  desired  pressure  profiles  were 
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scaled  to  prevent  saturation  of  either  the  source  array  D/A  or  the  feedback  array 
A/D.  Both  the  Kalman  filter  and  CNLLS  estimators  were  used  in  these  experiments. 
The  Kalman  filter  used  a  =  1.0  and  Pf  =  10“®I;  the  CNLLS  used  a  forgetting  factor 
of  7  =  0.9550  and  limited  the  condition  number  of  to  200.  This  limit  on 

the  condition  number  is  much  more  severe  than  the  10®  used  in  the  simulations  of 
Chapter  3  due  to  numerical  sensitivities  in  the  eigenvalue  library  implemented  in 
single  precision  on  the  TMS320C40.  Experiments  described  later  in  this  section  ex¬ 
plore  other  parameter  choices  for  the  channel  estimators  and  confirm  that  although 
these  may  not  be  the  ideal  parameter  choices  for  the  estimators  they  are  reasonable. 
Figure  4-4  shows  the  performance  of  the  feedback  control  algorithm  exciting  system 
mode  one.  Both  estimators  give  a  very  high  SER:  the  Kalman  filter  estimator  attains 
nearly  50  dB  of  SER,  while  the  CNLLS  estimator  levels  out  around  45  dB  SER.  The 
Kalman  filter  experiment  shows  several  drop-outs  of  the  signal.  These  were  caused 
by  network  delays  during  the  experiment.  The  transmitter  and  receiver  modems 
used  a  semaphore  protocol  via  files  on  the  NFS  server  to  synchronize  the  sampling 
of  the  pressure  field  at  the  feedback  array.  Unfortunately,  network  delays  sometimes 
disrupted  this  synchronization  and  the  receiver  did  not  sample  the  pressure  until 
after  the  source  array  stopped  transmitting.  These  glitches  are  regrettable,  but  they 
do  allow  the  Kalman  filter  estimator  to  demonstrate  its  robustness  to  this  sort  of 
error.  Other  experiments  shown  later  in  the  chapter  indicate  the  CNLLS  and  LMS 
algorithms  are  also  robust  to  this  sort  of  error.  This  is  encouraging  for  the  ocean 
experiments  because  the  radio  telemetry  link  shown  in  Fig.  1-1  between  the  feedback 
array  and  source  ship  will  likely  experience  similar  intermittent  disruptions. 

Figures  4-5  through  4-7  show  the  performance  of  both  estimators  for  system 
modes  two  through  four.  All  of  these  experiments  demonstrate  that  the  feedback 
algorithm  is  capable  of  exciting  high  fidelity  modes.  The  fidelity  of  the  excited 
pressure  profile  for  mode  four  is  particularly  impressive,  since  the  singular  value 
corresponding  to  this  system  mode  is  more  than  20  dB  below  the  singular  value 
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Figure  4-4;  SER  Performance  for  System  Mode  1  in  Time-Invariant  Waveguide 

for  mode  one.  This  can  be  seen  in  Fig.  4-2,  where  the  magnitude  of  the  pressure 
profile  of  system  mode  four  is  far  below  that  of  system  mode  one.  The  performance 
for  mode  four  shown  in  Fig.  4-7  indicates  this  mode  is  both  numerically  reachable 
and  observable  for  the  conditions  in  the  flume.  Figure  4-8  shows  the  performance  of 
the  feedback  control  algorithm  when  the  desired  pressure  profile  was  system  mode 
one  plus  0.4j  times  system  mode  two.  The  high  SER  for  this  experiment  confirms 
the  algorithm  can  excite  linear  combinations  of  the  system  modes  as  well  as  the 
individual  modes.  If  the  desired  acoustic  mode  in  an  ocean  experiment  is  not  one 
of  the  system  modes  but  falls  within  the  appropriate  numerical  span  of  the  system 
modes,  the  feedback  control  algorithm  should  be  able  to  excite  the  pressure  profile 
of  the  acoustic  mode  with  high  fidelity. 

Figures  4-9  through  4-11  show  the  performance  of  the  channel  estimators  with 
different  choices  of  parameters.  All  the  experiments  shown  in  these  figures  started 
with  no  a  priori  knowledge  of  Q,  i.e.,  random  Q[0],  and  attempted  to  excite  system 


103 


SER  (dB) 


45 


Figure  4-5:  SER  Performance  for  System  Mode  2  in  Time-Invariant  Waveguide 


Figure  4-6:  SER  Performance  for  System  Mode  3  in  Time-Invariant  Waveguide 
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Figure  4-7:  SER  Performance  for  System  Mode  4  in  Time-Invariant  Waveguide 


Figure  4-8:  SER  Performance  for  Mode  1  -f-  QAj  Mode  2  in  Time-Invariant  Waveg¬ 
uide 
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mode  one.  Figure  4-9  shows  the  performance  for  the  Kalman  filter  estimator  for 
various  values  of  Pf  when  a  =  1.0.  The  mean  magnitude  of  the  elements  of  Q 
is  0.16  for  this  experiment.  All  of  the  covariance  matrices  for  the  experiments  in 
Fig.  4-9  are  far  below  the  steady  value  of  most  elements  of  the  replica  matrix. 

With  a  =  1.0,  the  covariance  of  the  Kalman  filter  estimator  will  grow  in  all  direc¬ 
tions  except  w.  This  may  explain  why  the  smallest  Pf  gives  the  best  performance, 
as  the  estimator’s  covariance  matrix  will  grow  the  slowest  for  this  small  innovation 
covariance.  From  examining  this  figure,  it  is  tempting  to  suggest  the  high  frequency 
content  of  the  SER  increases  as  Pf  decreases.  Plotting  the  spectra  of  these  curves 
reveals  that  the  Pf  =  2  x  10~®I  trial  has  more  high  frequency  content  in  the  SER 
than  the  Pf  =  10“®I  trial.  More  importantly,  SER  is  not  necessarily  a  good  indicator 
of  the  frequency  content  of  the  errors.  Consider  the  case  when  the  observed  pressure 
p[n]  =  Pd  +  Ape-^‘^".  This  would  have  constant  SER  plotted  against  time  regardless 
of  the  frequency  u)  in  the  error  term.  Thus,  the  frequency  content  of  the  pressure 
error  may  not  be  refiected  in  the  SER.  It  is  also  important  to  remember  that  the 
curves  in  Fig.  4-9  are  from  a  single  trial.  The  limited  time  available  for  experiments 
precluded  running  extensive  trials  in  the  fiume  to  allow  averaging  of  performance 
over  many  trials. 

Figure  4-10  plots  the  performance  of  the  CNLLS  estimator  with  different  choices 
for  the  forgetting  factor  7.  The  relatively  smallest  forgetting  factor  (7  =  0.9,  effec¬ 
tive  window  length  =  22  samples)  converges  the  fastest  of  the  experiments  shown, 
although  7  =  0.925  reaches  a  slightly  higher  final  SER.  The  strong  performance  of 
these  relatively  short  windows  is  not  surprising  because  the  noise  level  in  the  flume 
is  so  low.  Consequently,  the  channel  estimator  does  not  need  to  average  many  it¬ 
erations  to  minimize  the  effect  of  the  observation  noise.  The  tighter  limit  on  the 
inverse  of  4^[n]  introduces  more  bias  into  the  channel  estimate  than  was  seen  in  the 
simulation  results  in  Chapter  3.  This  explains  why  the  CNLLS  performance  does 
not  equal  the  SER  obtained  by  the  Kalman  filter  in  the  waveguide,  although  the 
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CNLLS  often  performed  slightly  better  than  the  Kalman  filter  in  the  simulations  of 
time-invariant  ocean  channels. 


Figure  4-9:  Comparison  of  the  SER  for  System  Mode  1  in  the  Time-Invariant  Waveg¬ 
uide  for  different  Pf  in  the  Kalman  Filter  Channel  Estimator 

Figure  4-11  shows  the  SER  of  the  LMS  channel  estimator  for  two  different  choices 
of  a.  The  algorithm  converges  faster  and  to  a  higher  SER  for  a  =  0.5,  although  it 
does  not  match  the  performance  of  the  Kalman  filter  or  CNLLS  estimators.  It  is 
reassuring  that  these  latter  estimators  do  obtain  better  performance  than  the  LMS 
for  the  additional  computational  complexity  they  require.  Figure  4-12  plots  the 
pressure  magnitude  profile  against  depth  for  all  three  algorithms,  as  well  as  the 
desired  profile  for  mode  one.  The  Kalman  filter  and  CNLLS  estimators  generate  a 
pressure  profile  which  is  almost  indistinguishable  from  the  desired  pressure  profile. 
The  LMS  estimator  does  not  match  the  desired  profile  quite  as  closely,  but  does  still 
give  a  good  approximation  to  the  profile  of  system  mode  one. 
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Figure  4-10:  Comparison  of  the  SER  for  System  Mode  1  in  the  Time-Invariant 
Waveguide  for  different  7  in  the  CNLLS  Channel  Estimator 


Figure  4-11:  Comparison  of  the  SER  for  System  Mode  1  in  the  Time-Invariant 
Waveguide  for  different  choices  of  a  in  the  LMS  Channel  Estimator 
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Figure  4-12:  Depth  Profile  of  Pressure  Magnitude  for  All  Channel  Estimators  after 
100  Iterations 

4.4  Time- Varying  Channels 

The  experiments  presented  in  this  section  investigate  the  transient  performance  of  the 
feedback  control  algorithm.  These  experiments  consisted  of  allowing  the  algorithm 
to  iterate  for  a  fixed  amount  of  time,  then  inserting  a  scatterer  to  change  the  acoustic 
propagation  in  the  waveguide.  The  scatterer  used  was  a  40  x  19.5  x  19.5  cm  cinder 
block.  The  block  was  suspended  in  a  manner  allowing  it  to  be  consistently  placed 
at  r  =  3.15  m  in  the  waveguide  such  that  the  bottom  of  the  block  was  at  a  depth 
of  roughly  z  =  25  cm.  The  dynamic  characteristics  of  the  insertion  were  hand- 
controlled  by  a  human  operator,  and  were  not  as  consistent  between  trials  as  the 
location  was.  This  difficulty  resulted  in  the  two  strategies  for  transient  experiments 
outlined  in  the  following  sections. 

All  the  experiments  attempted  to  excite  system  mode  one  starting  from  a  cor¬ 
rupted  version  of  the  true  Q.  The  Kalman  filter  estimator  used  Pf  =  lO"^!  and 
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Oi  —  1.0,  while  the  CNLLS  estimator  had  7  =  0.95  (effective  window  length  ft;  45 
iterations)  and  a  condition  number  limit  of  200.  The  normalized  LMS  estimator 
used  a  =  0.5.  The  physical  configuration  of  the  waveguide  and  arrays  were  the  same 
as  described  earlier  in  this  chapter. 


4.4.1  Abrupt  Channel  Change 

These  experiments  changed  the  acoustic  properties  of  the  waveguide  such  that  from 
the  estimators’  point  of  view  the  change  occurred  instantaneously.  To  implement 
this,  the  feedback  control  algorithm  ran  for  100  iterations  with  an  open  waveguide, 
then  paused.  The  cinder  block  was  lowered  into  the  waveguide,  the  water  waves 
excited  in  the  tank  were  allowed  to  calm  down,  and  then  the  algorithm  resumed 
with  the  same  state  and  desired  pressure  profile  it  had  when  it  paused.  From  the 
estimators’  vantage,  the  brick  appeared  instantaneously  between  iterations  without 
disturbing  the  water.  This  experiment  was  repeated  for  each  of  the  channel  esti¬ 
mators.  This  method  of  modifying  the  acoustic  properties  of  the  channel  allowed 
the  experiments  to  focus  on  the  transient  behavior  of  the  estimators  independently 
of  the  settling  time  of  any  wave  dynamics  introduced  by  the  lowering  of  the  cinder 
block. 

Figure  4-13  presents  the  SER  curves  for  these  experiments.  The  presence  of  the 
block  does  alter  the  replica  matrix  and  thus  the  modes  of  the  system.  The  desired 
pressure  profile  at  the  feedback  array  is  kept  unchanged  as  the  first  mode  of  the  open 
waveguide.  Projecting  this  pressure  profile  on  the  first  three  modes  of  the  replica 
matrix  with  the  block  in  place  reveals  more  than  99%  of  the  energy  of  the  desired 
profile  still  falls  in  the  span  of  the  first  three  system  modes,  i.e.. 
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where  uq^j  is  the  z*'*’  column  of  the  replica  matrix  Q  with  the  block  inserted  in 
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the  waveguide,  and  pj.  is  the  first  mode  of  the  system  for  the  open  waveguide.  The 
dotted  line  at  roughly  30  dB  indicates  the  best  possible  performance  exciting  mode 
one  with  the  block  in  place.  This  bound  on  the  SER  is  down  from  a  theoretical 
limit  of  150  dB  in  the  absence  of  noise  for  the  open  tank  without  the  block.  As  the 
observation  noise  limits  the  performance  to  45-50  dB  of  SER,  inserting  the  block 
changes  the  environment  from  one  where  the  performance  is  noise-limited  to  one 
where  the  performance  is  propagation  limited  by  the  coupling  induced  by  the  block. 
The  Kalman  filter  estimator  adapts  to  the  new  channel  in  roughly  12  iterations. 
The  CNLLS  and  LMS  estimators  take  longer  to  converge  after  the  change  in  the 
channel.  The  CNLLS  estimator  takes  about  40  iterations  to  converge,  while  the 
LMS  algorithm  hasn’t  converged  100  iterations  after  the  block  appears.  Both  of  these 
convergence  characteristics  could  probably  be  improved  by  modifying  the  estimator 
parameters:  decreasing  7  for  the  CNLLS  or  increasing  jl  for  the  normalized  LMS. 
While  these  changes  would  increase  the  convergence  rate,  they  would  also  make  the 
estimators  more  sensitive  to  observation  noise.  Decreasing  7  would  decrease  the 
effective  window  length  the  CNLLS  estimator  averages  over,  while  a  larger  jl  makes 
the  LMS  estimator  more  sensitive  to  noise  even  after  it  has  converged.  In  noisy 
ocean  experiments,  this  sensitivity  would  be  a  more  serious  issue  than  in  the  very 
quiet  and  controlled  laboratory  waveguide  setting. 

4.4.2  Gradual  Channel  Change 

In  the  experiments  described  in  this  section,  the  feedback  control  algorithm  ran  con¬ 
tinuously  while  the  cinder  block  was  lowered  into  the  waveguide,  rather  than  pausing 
as  in  the  experiments  described  in  the  previous  section.  The  planned  sequence  was 
to  allow  the  algorithm  to  run  for  50  iterations,  then  gradually  insert  the  block  over 
the  next  ten  to  fifteen  iterations  while  trying  to  minimize  the  water  waves  excited 
in  the  tank  by  the  disturbance.  The  execution  was  not  always  successful  in  terms  of 
the  timing  of  the  insertion  or  minimizing  the  fluid  dynamic  transient. 
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Figure  4-13:  Transient  Response  of  Channel  Estimators  for  Abrupt  Change  in  Chan¬ 
nel 

Figures  4-14  through  4-16  present  the  results  of  these  experiments.  The  different 
line  types  on  each  figure  represent  successive  trials  of  the  same  experiment.  It  is  not 
easy  to  extract  the  sort  of  quantitative  results  obtained  in  the  previous  section  from 
these  plots.  The  Kalman  filter  clearly  converges  the  fastest  in  these  experiments, 
with  the  CNLLS  and  LMS  estimators  following  in  that  order.  Both  the  Kalman 
filter  and  CNLLS  estimators  come  very  close  to  the  LS  bound  on  SER  with  the 
block  in  the  waveguide.  This  bound  for  these  trials  is  not  necessarily  the  same  as 
the  one  shown  in  Fig.  4-13  since  the  waveguide  characteristics  drift  and  these  trials 
were  performed  at  a  different  time  and  took  longer  to  complete.  Consequently,  the 
extent  to  which  Pd  fell  within  the  numerical  span  of  Q  varied  from  trial  to  trial.  This 
bound  only  varied  by  a  few  dB  over  the  experiments  shown  in  this  section. 

The  longer  convergence  times  of  all  the  estimators  in  these  experiments  when 
compared  to  Fig  4-13  could  be  due  to  a  combination  of  several  factors.  The  acoustic 
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properties  of  the  waveguide  are  changing  every  iteration  for  roughly  a  dozen  iterations 
due  to  the  lowering  of  the  block  into  the  waveguide.  In  addition,  the  water  waves 
excited  by  inserting  the  block  will  also  affect  the  acoustic  propagation.  Even  a  very 
gentle  touch  lowering  the  block  into  the  channel  excited  water  waves  that  took  several 
iterations  to  settle. 


Figure  4-14:  Transient  Response  of  Kalman  Filter  Channel  Estimator  for  Gradual 
Change  in  the  Waveguide 

A  rough  quantitative  estimate  of  the  increase  in  convergence  time  can  be  obtained 
by  starting  from  the  observation  that  in  most  of  the  trials  in  these  experiments  the 
block  was  completely  lowered  by  iteration  65.  Looking  closely  at  Fig.  4-14  reveals  the 
Kalman  filter  estimator  generally  converged  within  25  iterations  after  the  block  was 
in  place.  This  was  about  twice  as  long  as  the  convergence  time  observed  in  Fig.  4-13. 
A  similar  examination  of  Fig.  4-15  reveals  an  increase  of  convergence  time  to  roughly 
70  iterations,  rather  than  the  40  iterations  observed  in  Fig.  4-13.  The  LMS  algorithm 
never  completely  converged  in  either  set  of  experiments,  making  it  impossible  to 
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Figure  4-15:  Transient  Response  of  CNLLS  Channel  Estimator  for  Gradual  Change 
in  the  Waveguide 


Figure  4-16:  Transient  Response  of  LMS  Channel  Estimator  for  Gradual  Change  in 
the  Waveguide 
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make  a  similar  comparison  for  this  estimator.  Comparing  these  transient  response 
times  to  those  seen  with  the  abrupt  channel  change,  one  qualitative  conclusion  that 
can  be  made  is  water  dynamics  slow  the  convergence  of  the  estimators.  While  the 
insertion  of  the  block  was  not  perfectly  consistent  for  each  trial,  the  water  waves 
were  excited  in  roughly  the  same  way.  The  varying  transient  response  times  among 
the  estimators  indicate  that  the  estimators  experience  differing  levels  of  difficulty 
tracking  the  water  dynamics.  If  the  dynamics  were  slow  enough  that  the  estimators 
could  perfectly  track  the  changes  in  Q,  all  of  the  transient  times  of  the  estimators 
would  be  the  settling  time  of  the  water  dynamics.  The  variation  of  the  response 
times  in  Figs.  4-14  through  4-16  confirm  the  estimators  are  unable  to  keep  up  with 
the  changes  in  Q  due  to  the  hydrodynamics.  The  Kalman  filter  does  not  assume 
the  replica  matrix  dynamics  are  a  stationary  process,  so  it  is  not  surprising  that 
it  adapts  most  quickly.  In  fact,  the  transient  time  for  the  Kalman  filter  might  be 
taken  as  a  rough  estimate  of  the  settling  time  for  the  dynamics  of  the  water  waves. 
The  CNLLS  and  LMS  estimators  assume  Q  is  stationary,  so  they  won’t  be  able  to 
make  substantial  progress  converging  on  the  new  Q  with  the  block  in  place  until  the 
non-stationarity  introduced  by  the  water  waves  dies  out.  For  the  CNLLS  estimator, 
the  increase  in  convergence  time  from  the  abrupt  channel  change  experiment  should 
reflect  the  settling  time  for  the  water  dynamics.  The  increase  in  convergence  time 
of  roughly  30  iterations  is  commensurate  with  the  estimate  of  the  time  constant  for 
the  hydrodynamics  obtained  from  the  Kalman  filter  transient.  The  duration  of  the 
ripples  in  the  LMS  SER  in  Fig.  4-16  also  supports  this  estimate  of  the  hydrodynamics 
time  constant.  As  noted  in  the  discussion  of  Fig  4-9,  time  characteristics  of  the  SER 
do  not  necessarily  reveal  all  the  underlying  dynamics  of  the  error,  but  here  the  SER 
curve  for  the  LMS  estimator  does  corroborate  the  estimated  settling  time  for  the 
water  waves  obtained  from  the  CNLLS  and  Kalman  filter  SER  curves. 

The  estimators’  difficulty  tracking  the  flume  water  dynamics  are  not  necessarily 
cause  for  concern  in  the  ocean  environment.  Visual  observation  indicated  the  waves 
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excited  by  lowering  the  block  propagated  the  8  m  from  the  block  to  the  feedback  array 
in  less  than  3  seconds.  Scaling  up  by  a  factor  of  20  for  a  400  Hz  ocean  experiment,  this 
indicates  that  in  an  ocean  experiment  the  estimators  could  not  track  water  dynamics 
propagating  26  m/s,  an  excessive  velocity  for  ocean  phenomena.  The  simulations  of 
the  solitary  wave  in  Section  3.4  suggests  the  Kalman  filter  can  track  ocean  processes 
propagating  at  least  as  fast  as  1  m/s,  which  is  a  generous  estimate  for  the  group 
velocity  of  these  waves. 

4.5  Summary 

This  chapter  presented  a  sequence  of  experiments  using  the  proposed  feedback  con¬ 
trol  algorithm  in  a  laboratory  waveguide.  These  experiments  demonstrate  that  the 
algorithm  performs  well  in  time-invariant  channels  when  implemented  on  real-time 
signal  processing  hardware  connected  to  real  sources  and  hydrophones.  This  en¬ 
vironment  provided  a  very  useful  test  bed  for  verifying  the  algorithm  was  feasible 
and  making  some  preliminary  investigations  of  its  transient  behavior.  More  de¬ 
tailed  investigations  of  the  tracking  performance  of  the  estimators  would  require  a 
well-controlled  means  of  introducing  a  fluid  dynamic  disturbance  of  the  appropriate 
bandwidth  in  the  waveguide.  The  results  obtained  in  the  experiments  described  in 
this  chapter  are  encouraging  for  the  success  of  an  ocean  experiment.  In  addition,  the 
experiments  provided  some  valuable  experience  for  the  deployment  and  debugging 
of  a  real-time  system  implementing  the  feedback  control  algorithm. 
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Chapter  5 


Conclusions  and  Future  Directions 


The  work  presented  in  this  thesis  develops  an  adaptive  feedback  controller  to  excite 
a  single  mode  in  the  shallow  water  channel.  No  mid-frequency  (circa  400  Hz)  ocean 
experiment  has  successfully  generated  a  single  mode  in  this  environment.  A  single 
mode  source  in  this  frequency  range  would  be  very  useful  for  studying  oceanographic 
processes  in  the  coastal  environment.  The  control  algorithm  incorporates  elements  of 
control  theory,  adaptive  estimation,  array  processing  and  underwater  acoustics.  The 
controller  iterates  between  choosing  the  source  array  weight  vector,  and  updating 
its  estimate  of  the  replica  matrix  between  the  source  and  feedback  arrays.  Three 
channel  estimators  are  developed  based  on  different  models  for  the  properties  and 
statistics  of  the  replica  matrix. 

The  work  leading  to  the  control  algorithm  also  yielded  two  additional  results. 
The  first  is  a  formulation  of  the  maximum  a  posteriori  (MAP)  mode  filter.  This 
filter  segues  between  the  pseudo-inverse  and  sampled  mode  shape  filters  as  the  hy¬ 
drophone  array  aperture  decreases.  This  result  provides  a  theoretical  basis  for  an 
earlier  ad  hoc  “dropped  eigenvalue”  mode  filter  proposed  by  Yang  [41].  The  MAP 
mode  filter  also  has  an  appealing  interpretation  as  a  generalization  of  the  discrete 
spatial  Wiener  filter.  The  investigation  of  controllability  and  observability  provides 
a  cohesive  framework  for  examining  numerical  issues  germane  to  a  practical  exper- 
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iment.  Specifically,  criteria  are  given  for  determining  when  a  desired  mode  can  be 
excited  and  observed  given  proposed  array  configurations,  observation  noise  charac¬ 
teristics,  source  power  limitations  and  approximate  acoustic  attenuation  in  range. 

A  series  of  simulations  presented  in  Chapter  3  investigate  the  performance  of  the 
proposed  algorithm  in  several  realistic  shallow  water  environments:  a  range-invariant 
waveguide,  a  rock  outcropping,  a  downsloping  wedge,  and  a  propagating  solitary  in¬ 
ternal  wave.  The  performance  of  the  various  channel  estimators  in  each  environment 
are  compared  to  each  other  and  the  performance  of  an  open  loop  controller.  The 
Kalman  filter  and  CNLLS  estimators  both  generated  high  fidelity  modes  in  these 
simulations,  with  the  latter  usually  performing  slightly  better. 

The  proposed  algorithm  was  also  evaluated  in  a  series  of  laboratory  waveguide 
experiments.  The  algorithm  successfully  excited  several  different  modes  of  the  sys¬ 
tem.  The  Kalman  filter  channel  estimator  proved  more  robust  to  the  computational 
limitations  of  the  real-time  DSP  hardware  in  these  experiments.  The  waveguide 
experiments  also  provided  some  preliminary  data  on  the  transient  response  of  the 
channel  estimators.  The  Kalman  filter  also  performed  better  than  the  other  estima¬ 
tors  in  the  transient  experiments.  The  success  of  the  algorithm  in  the  laboratory 
waveguide  is  encouraging  for  the  prospects  of  an  ocean  experiment. 

Several  aspects  of  the  single  mode  excitation  problem  present  intriguing  directions 
for  future  research.  One  possible  improvement  in  the  feedback  control  algorithm 
would  be  in  the  weight  vector  selection  step.  Rather  than  choosing  the  source  array 
weight  vector  w[n]  to  minimize  the  squared  error  between  the  predicted  pressure 
Q[n]w[n]  and  desired  pressure  vector  p^,  w[n]  could  be  chosen  to  minimize  the  ratio 
of  the  energy  in  all  modes  to  the  energy  in  the  desired  mode.  In  the  limit  with  only 
mode  mo  excited,  this  ratio  would  be  1.  More  precisely,  the  optimization  problem 
the  new  weight  selection  step  would  solve  is 


w[n 


arg  min 


w^B^[n]B[n]w 

[n]w’ 


(5.1) 
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where  B[n]  is  an  estimate  of  the  source  weight  to  mode  coefficient  transfer  matrix 
from  Eq.  1.35  and  bj„oM  is  the  row  of  that  estimate.  As  formulated  in  Eq.  5.1, 
the  optimization  does  not  have  a  unique  solution,  but  by  introducing  an  additional 
constraint,  a  unique  solution  can  be  obtained.  The  constrained  optimization  is 

w[n]  =  argimnw^  B«[n]  B[n]w  subject  to  b^oMw  =  1.  (5.2) 

and  it  is  clear  any  solution  to  this  optimization  will  also  be  a  solution  to  Eq.  5.1. 
The  constrained  problem  in  Eq.  5.2  has  the  same  form  as  the  Minimum  Variance 
Spectral  Estimator  [66],  [40].  The  solution  to  this  optimization  can  be  written  in 
closed  form; 

(BH[n]B[n]) 

w[n]  =  \  - :  .-1. - •  (5.3) 

bmoN  (B«‘[n]B[n]) 

The  performance  of  this  weight  selection  criteria  would  need  to  be  evaluated  in  simu¬ 
lations  and  laboratory  waveguide  experiments  and  compared  against  the  performance 
of  the  least-squares  criteria  used  in  the  thesis. 

The  proposed  algorithm  uses  a  very  basic  model  for  the  acoustic  propagation 
in  the  feedback  volume  and  also  a  simple  feedback  control  algorithm.  The  simu¬ 
lations  and  laboratory  experiments  indicate  these  straightforward  approaches  can 
successfully  excite  a  single  mode.  This  success  encourages  further  investigations  in¬ 
corporating  more  sophisticated  control  algorithms  and  ocean  models  to  see  if  these 
advances  can  provide  even  better  performance.  Specifically,  the  field  of  robust  con¬ 
trol  contains  many  important  results  that  may  provide  additional  insight  into  the 
single  mode  control  problem  [49],  [50],  [51].  In  pursuing  these  advances,  a  balance 
must  be  made  between  improved  performance  exciting  a  mode  and  greater  sensi¬ 
tivity  to  mismatch  errors  in  either  the  ocean  acoustics  or  estimated  plant.  While 
the  algorithm  proposed  here  is  very  modest  in  its  technical  scope,  it  is  flexible  in 
responding  to  a  wide  range  of  ocean  environments.  Any  more  sophisticated  control 
algorithm  will  want  to  preserve  this  robustness  as  much  as  possible. 
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Another  possible  improvement  would  be  a  more  aggressive  temporal  sampling 
strategy  for  the  feedback  volume.  The  iteration  time  could  be  reduced  from  the 
propagation  time  of  the  slowest  mode  between  the  source  and  feedback  arrays  to  the 
difference  between  the  travel  times  of  the  fastest  and  slowest  modes.  This  would 
require  accurate  estimates  of  the  mode  group  velocities.  These  estimates  could  be 
obtained  numerically  from  the  sound  speed  profiles  at  the  source  and  feedback  ar¬ 
rays  [32]  and  if  the  estimates  are  sufficiently  accurate  and  reliable,  several  packets 
of  modes  could  be  en  route  between  the  source  and  feedback  arrays  simultaneously. 
This  decrease  in  the  iteration  time  for  the  control  algorithm  would  allow  the  channel 
estimator  to  track  a  higher  bandwidth  of  oceanographic  processes.  For  the  ranges 
and  environments  examined  in  this  thesis,  the  conservative  iteration  interval  set  by 
the  propagation  time  of  the  slowest  mode  sufficiently  samples  the  most  important 
processes,  so  the  incremental  gains  obtained  by  decreasing  the  iteration  time  do  not 
merit  the  significantly  increased  risk  of  poor  performance  if  the  estimated  group  ve¬ 
locities  are  incorrect.  Inaccurate  group  velocity  estimates  could  cause  mode  packets 
that  are  no  longer  temporally  distinct  at  the  feedback  array.  Future  deployments  of 
the  algorithm  may  occur  in  environments  with  ocean  processes  operating  on  shorter 
time  scales,  or  with  a  greater  distance  between  the  source  and  feedback  arrays.  In 
such  scenarios,  the  benefits  of  reducing  the  iteration  time  might  merit  the  additional 
complexity  and  risk  associated  with  this  strategy.  A  specific  instance  where  this 
could  prove  valuable  would  be  scenarios  where  surface  wave  processes  have  a  signif¬ 
icant  effect  on  the  mode  propagation.  Surface  waves  have  much  shorter  time  scales 
than  most  internal  ocean  processes  and  reducing  the  iteration  time  of  the  control 
algorithm  could  prove  crucial  to  allowing  the  estimators  to  track  these  processes. 

The  single  mode  source  provided  by  the  feedback  control  algorithm  can  be  used 
for  other  acoustic  measurements  besides  estimating  the  horizontal  wavenumber  spec¬ 
tra  of  an  internal  wave  packet.  When  the  propagation  in  the  observation  volume  is 
well-modeled  by  a  discrete  set  of  propagating  modes,  the  mode  coupling  over  that 
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volume  gives  a  simple  parametric  model  for  predicting  the  propagation  through  that 
volume.  Long-term  measurements  of  these  coupling  statistics  could  provide  use¬ 
ful  measurements  of  the  dynamics  of  acoustic  propagation  in  coastal  regions  which 
could  be  exploited  by  more  sophisticated  ocean  acoustics  signal  processing  algo¬ 
rithms.  Currently,  very  little  data  is  available  about  the  statistics  and  dynamic 
behavior  of  propagation  in  these  regions. 

The  feedback  control  algorithm  may  have  other  oceanographic  applications  be¬ 
sides  a  single  mode  source.  Gingras  [5]  proposed  a  single  mode  source  as  a  method 
of  minimizing  environmental  backscatter  for  an  active  sonar  system.  While  a  sin¬ 
gle  propagating  mode  is  a  useful  conceptual  approach  to  visualize  which  portions 
of  the  water  column  would  be  illuminated  by  such  a  system,  such  a  pressure  pro¬ 
file  may  not  be  the  optimal  transmitted  pressure  field  to  minimize  backscatter  from 
the  bathymetry.  The  feedback  control  algorithm  could  adapt  the  source  array  to 
minimize  the  backscattered  energy  in  the  absence  of  targets. 

Another  possible  but  highly  speculative  use  for  the  feedback  control  algorithm 
would  be  in  a  communication  scenario.  Different  messages  could  be  modulated 
onto  distinct  modes  for  a  communication  network.  This  would  allow  simultaneous 
communication  at  the  same  frequency  by  exploiting  the  modal  diversity  of  the  en¬ 
vironment.  Any  communication  system  proposed  utilizing  this  modulation  would 
have  to  address  many  issues  such  as  mode  coupling  due  to  range-inhomogeneities  in 
the  water  column  or  bottom  bathymetry  between  the  feedback  array  and  receiver. 
The  concept  is  attractive  and  bears  closer  examination. 

The  ideas  presented  in  this  section  are  intriguing  as  extensions  or  improvements 
on  the  algorithm  presented  in  the  thesis.  At  the  current  time,  the  most  pressing 
future  work  to  be  done  is  an  actual  ocean  test  of  the  algorithm  proposed.  Significant 
engineering  challenges  such  as  the  radio  telemetry  link  from  the  feedback  array  to 
the  source  ship  and  continuous  mode  shape  estimation  at  the  feedback  array  location 
must  still  be  addressed.  The  simulations  and  laboratory  experiments  indicate  that 
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the  feedback  control  algorithm  should  be  capable  of  exciting  a  single  mode.  If  this 
experiment  is  successful,  more  advanced  uses  can  be  explored.  If  the  trial  reveals 
unexpected  shortcomings  or  difficulties  with  the  feedback  control  algorithm,  these 
problems  can  be  addressed.  Without  such  an  experiment,  any  future  work  performed 
with  the  algorithm  must  be  viewed  as  speculative  at  best. 
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Appendix  A 


Transient  Response  of  the  CNLLS 
Estimator 


Section  2.3  describes  the  CNLLS  estimator  and  notes  that  it  is  important  to  limit 
the  condition  number  of  at  each  iteration  and  not  the  condition  number  of  $ 
propagated  by  the  algorithm.  Specifically,  the  estimator  may  respond  more  slowly 
to  abrupt  changes  in  the  desired  mode  or  the  channel  if  the  condition  number  of 
#  is  limited.  To  understand  why,  consider  the  simplified  case  where  w[n]  =  wq 
for  n  =  1,...,277,  where  t]  is  the  limit  on  the  condition  number.  The  eigenvalue 
decomposition  of  ^[2r]]  can  be  written 


where 


^[2rj]  = 


2ri  0  •••  0 

0  O'-.: 

:  0 

0  •••  0  0 


■ 

Wo 

V2 

Vl 

(A.1) 


(A.2) 
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Suppose  the  desired  mode  changes  abruptly  at  n  =  277  such  that  w[n]  is  now  or¬ 
thogonal  to  Wo  for  n  >  27].  Note  that  if  w[n]  is  orthogonal  to  Wq,  we  may  assume 
w[ra]  =  V2  without  loss  of  generality,  since  we  can  pick  the  Vj  to  be  any  orthonormal 
basis  of  orthogonal  complement  of  wo.  Moreover,  assume  QN  has  errors  that  were 
orthogonal  to  wq  but  are  not  orthogonal  to  V2,  i.e.,  (QW  -QN)v2  #  0.  The  change 
of  w[n]  to  V2  will  cause  Q[n]  to  converge  to  a  more  accurate  estimate  of  Q[n]  such 
that  (Q[n]  —  Q[n])v2  =  0. 

The  rate  of  this  convergence  depends  on  which  of  #  or  has  its  condition 
number  limited.  From  Eqs.  2.21  and  2.22,  the  correction  to  Q[n]  is  a[n]k^[n],  where 
a[n]  is  (p[n]  -  Q[n]w[n])  =  (Q[n]  -  Q[n])w[n]  and  k[n]  =  $“^[n]w[n]. 

Let  ki[n]  be  the  gain  vector  when  $[n]  is  propagated  with  its  condition  number 
limited,  and  let  k2[n]  be  the  gain  vector  when  $[n]  propagates  without  any  limit 
on  its  condition  number,  but  the  condition  number  is  limited  at  each  iteration  when 
^CL  is  computed  per  Eq.  2.25.  The  notation  {■)cl  denotes  the  operation  of  limiting 
the  condition  number  of  its  argument  to  be  77  while  leaving  the  first  eigenvalue 
unchanged.  Then 


ki[277  +  l]  =  [277  +  1])  w[277  +  1] 

T  -1 

277  0  •••  0 


V^^V2 


(A.3) 


(A.4) 


[  0  0  2 
=  V2/3. 


(A.5) 
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Contrastingly,  for  the  CNLLS  estimator  given  in  Sec.  2.3, 


k2[2ry  +  l]  =  [277  +  1]w[2?7  +  1] 


277 

0  • 

••  0 

1 

0 

0 

0 

0 

...  0 

0 

y 

=  V2/2 


(A.6) 


(A.7) 


(A.8) 


The  gain  vector  k[n]  has  the  same  direction  in  each  case,  but  takes  a  larger  step  for 
the  case  where  the  condition  number  is  propagated  without  limiting,  then  limited 
after  inversion.  If  the  correlation  matrix  $  is  propagated  with  condition  number 
limiting,  the  gain  of  the  correction  is  reduced  even  the  first  step  after  w[n]  changes. 
In  fact,  if  k2 [277 +  2]  will  still  be  V2/2,  while  ki[277  +  2]  =  V2/4.  As  long  as  w[77]  =  V2, 
ki[n]  and  k2[77]  will  have  the  same  direction,  but  ||ki[n]||  <‘||k2[n]||,  so  k2[77]  will 
have  a  faster  transient  response. 

The  scenario  here  simplifies  the  case  in  several  aspects,  but  the  underlying  results 
are  still  valid:  propagating  $[77]  with  its  condition  number  unlimited  then  limiting 
the  condition  number  at  the  inversion  for  each  iteration  gives  better  transient  per¬ 
formance  than  propagating  a  #[77]  whose  condition  number  is  limited. 
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